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ABSTRACT 


In  this  work,  the  control  laws  of  a  bank-to*turn  mis* 
sile  using  and  optimal  estimator  in  the  terminal  guidance 
phase  were  designed,  and  the  effect  of  increasing  the  number 
of  measurement  sensors  in  the  missile  to  generate  more 
information  on  the  state  was  investigated.  In  the  design  of 
the  control  law  the  modern  optimal  control  theory  with  the  a 
quadratic  performance  index  was  used.  Implementation  of 
this  control  law  required  the  use  of  a  Kalman  filter  as  the 
optimal  estimator.  The  extended  Kalman  filter  algorithm  was 
utilized  in  the  present  study  since  the  measurement  states 
were  non-linear  functions  of  the  state  vectors.  In  order  to 
test  the  effects  of  the  implementation  of  the  increased  mea¬ 
surement  sensors,  two*,  four-,  and  six-measurement  sensors 
were  assumed  to  be  implemented  in  the  optimal  estimator.  By 
computer  analysis,  the  designed  guidance  laws  were  evaluated 
and  the  effect  of  the  implementation  of  increased  measurement 
sensors  was  tested. 

The  results  of  the  simulation  revealed  that  the  designed 
guidance  law  was  successful  within  the  specified  scenarios, 
the  effect  of  the  implementation  of  increased  measurement 
sensors  for  the  estimator  was  favorable  only  in  that 
increased  measurement  sensors  generated  more  information 
about  the  state  vectors. 
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I.  INTRODUCTION 


The  bank-to-turn  missile  has  high  lift  acceleration  in 
a  direction  perpendicular  to  its  wings.  For  airbreathing 
missiles  which  are  required  for  large  stand-off  ranges  it 
also  offers  the  advantage  of  lower  inlet  angles  of  attack 
than  that  of  skid-to-turn  missiles. 

Although  the  control  laws  using  the  modern  guidance 
control  theory  are  more  complex  than  the  normally  used 
proportional  navigation  air-to-ground  method,  they  have 
great  potential  for  maneuvering  targets  in  air-to-air 
engagement  situations.  The  application  of  optimal  control 
and  estimation  theory  to  bank-to-turn  missile  configurations 
has  been  tried  in  order  to  obtain  increased  performance.  In 
the  application  of  the  optimal  control  theory  to  the  bank- 
to-turn  missile,  it  is  necessary  to  have  information  on  all 
states  or  an  estimate  of  the  state  variables  that  have  not 
been  measured.  Since  the  state  information  available  from 
the  typical  missile  sensors  is  limited,  it  is  necessary  to 
employ  an  estimator.  The  estimator  used  in  this  study  was 
a  first  order  extended  Kalman  filter. 

The  present  work  addressed  the  design  and  evaluation  of 
the  optimal  state  estimator  and  optimal  control  laws  for 
application  to  a  bank-to-turn  missile.  In  the  development 
of  this  work,  Chapter  2  will  describe  the  kinematics  of  the 
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missile  and  formulate  the  equations  of  the  motion.  All 
major  assumptions  are  listed.  The  optimal  controller  is 
described  in  Chapter  3.  The  simulation  results  of  the 
controller  on  three  scenarios  described  later  are  used  to 
check  the  suitability  of  the  optimal  controller.  The 
optimal  estimator  is  covered  in  Chapter  4,  the  formulation 
of  the  state  equation,  measurement  equations  and  the 
derivations  of  the  elements  of  Jacobean  matrix  are  covered. 
The  extended  Kalman  filter  algorithm  for  non-linear  system 
is  then  reviewed.  Before  the  final  computer  simulation  a 
discussion  is  given  on  the  nominal  parameters  needed  in 
initializing  the  Kalman  filter.  The  results  of  the  simula¬ 
tion  of  estimators  with  two-7  four-,  and  six-measurement 
vectors  on  a  typical  scenario  are  then  presented.  In 
Chapter  5,  the  simulation  results  of  the  control  laws  imple¬ 
mented  the  estimator  with  two-,  four-  and  six-measurement 
sensors  on  three  scenarios  is  described.  The  detail  analy¬ 
sis  of  the  results  is  developed  in  order  to  investigate  the 
effect  of  the  key  variables  of  the  control  system  on  a  bank- 
to-turn  missile.  The  mean  miss  distance  determined  from  50 
Monti  Carlo  runs  as  a  performance  standard  was  used  in  the 
analysis  of  the  results  of  the  control  laws  as  a  function  of 
the  key  variables.  Finally,  the  conclusion  are  summarized 
in  the  last  Chapter. 


V. 


II.  KINEMATICS  OF  A  BANK-TO-TURN  MISSILE 


A.  ASSUMPTIONS 

The  geometry  and  the  equations  of  motion  of  the  missile 
will  be  developed,  representing  the  positions  in  space, 
under  the  following  assumptions: 

1.  The  local  geographic  coordinate  system  will  be  used 
as  the  inertial  reference. 

2.  To  simplify  the  equation  of  motion  the  stability 
axes  are  adopted  for  the  missile  body. 

3.  The  missile  velocity  due  to  thrust  is  constant. 
There  is  no  missile  acceleration  in  velocity  due  to 
thrust . 

4.  Each  control  surface  or  surfaces  is  rigid. 

5.  Relative  to  the  body  axes,  each  control  system  has 
only  one  degree  of  freedom.  The  missile's  control 
acceleration  vector  acts  normal  to  the  velocity 
vector,  that  is,  the  dot  product  of  acceleration 
vector  and  velocity  vector  is  zero,  also  the 
acceleration  acts  through  the  missile  center  of 
gravity . 

6.  The  Euler  angles  $  ,  6v  and  <t»  will  be  used  to 

describe  the  orientation  of  the  missile  with 
respect  to  inertial  space,  where  <*v  and  are 

horizontal  and  vertical  flight  path  angles,  and  $ 

is  the  roll  or  bank  angle. 

IS 


7.  The  first  order  lags  with  time  constants  x  and 

P 

tq  to  a  roil  rate  command  Pc  and  a  normal 
acceleration  command  ac  will  be  considered  for 
the  dynamic  response  of  the  missile. 

B.  MISSILE  CONTROL  METHOD 

Before  going  into  the  mathematical  detail  concerning 
the  motion  of  a  missile  in  space  as  a  result  of  a  guidance 
command,  it  is  helpful  to  review  with  the  control  law  for  a 
bank-to-turn  missile.  The  guidance  system  detects  whether 
the  missile  is  flying  too  much  or  too  high  to  the  left, 
right  or  vertical.  The  guidance  system  measures  these 
deviations  or  errors  and  sends  signals  to  the  control  system 
to  reduce  these  errors  to  aero.  The  task  of  the  control 
system  therefore  is  to  maneuver  the  missile  quickly  and 
efficiently  as  a  result  of  these  signals.  In  a  bank-to-turn 
missile  system,  the  guidance  angular  error  detector  produces 
two  signals  R  and  <t>  in  terms  of  polar  coordinate  expres¬ 
sion,  showed  in  Figure  2.1.  The  same  signals  can  be 
expressed  in  another  way,  that  is,  in  cartesian  coordinate. 
The  usual  method  is  to  regard  the  <p  signal  as  a  command 
to  roll  through  an  angle  a  measured  from  the  vertical 
and  then  to  maneuver  outwards  by  means  of  the  missile's 
elevators.  The  method  of  maneuvering  the  missile  is  as 
follows.  The  4>  command  goes  as  a  positive  command  to  one 
control  surface  and  a  negative  command  to  the  other,  this 


FIGURE  2.1.  MISSILE  CONTROL  SYSTEM  IN  POLAR  CO-ORD. 

causes  the  missile  to  roll.  The  R  command  goes  to 
surfaces  always  as  a  positive  demand,  this  causes  the 
missile  to  accelerate  normal  to  the  velocity  vector.  The 
intention  is  to  make  the  response  in  roll  fast  so  that  the 
commands  can  be  applied  simultaneously  which  makes  for 
simplicity,  only  a  pair  of  control  surfaces  are  used  as 
ailerons  and  elevators  at  the  same  time,  control  is  obtained 
by  means  of  the  separated  servos. 

C.  GEOMETRY  OF  MOTION  IN  SPACE 

The  motion  of  a  missile  as  a  particle  may  be  described 
by  using  coordinate  measured  with  moving  axes  (relative- 
motion  analysis) .  The  analysis  of  motion  can  be  simplified 
by  using  measurements  made  with  respect  to  a  moving  coordi¬ 
nate  system.  In  present  work,  the  equations  of  motion  will 
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be  described  with  respect  to  the  initial  missile  body  axes 
taken  as  the  inertial  axes.  After  launching,  the  missile 
body  axes  changes  position,  but  the  relative  position, 
velocity  and  acceleration  are  still  computed  in  inertial 
axes,  that  is,  in  the  initial  body  axes.  The  current  zero 
effort  miss  distance  is  calculated.  This  will  then  be 
transformed  to  values  in  missile  body  axes  to  compute  the 
amount  of  the  control  inputs  to  the  missile  using  Euler's 
transformation. 

Figure  2.2  shows  the  missile  and  target  as  points  with 
respective  vector  velocity  V  and  V't  .  In  this  analysis, 
because  |Vm|  is  assumed  to  be  at  least  2  j  V  |  ,  and  the 

angle  of  attack  is  assumed  to  be  small,  the  lead  angle  <p 
between  and  the  missile  longitudinal  axis  is  small. 

Let  the  vector  denote  the  relative  position  of  the  target 
with  respect  to  the  missile.  The  orientation  of  the  sight 
line  vector  in  inertial  space  can  be  represented  by  the 
angle  ^  and  0R  ,  as  shown  in  Figure  2.2.  The  sight 
line  vector  is  resolved  into  three  components  in  inertial 
axes  as  follows: 


Rx  3  R*cos (0^) *ccs  Opp)  =  Xt 

-  Xm 

(2.1) 

Ry  *  R*cos  (0R)  *sin(i{/R)  =  Yt 

-  Ym 

(2.2) 

Rz  *  -  R*sin  ($R)  =  Zt 

-  Zm 

(2.3) 

Where  R  is  the  magnitude  of  the  sight  line  vector. 
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Define  the  relative  velocity  vector: 


(2.4) 


Where  V*  and  V  are  target  and  missile  velocities, 
t  m  ° 

The  rate  of  change  of  the  sight  line  vector's  magnitude 
is  equal  to  the  component  of  V„  along  R  .  It  can  be 
expressed  as  following: 


R  = 


iyt  Rl 

R 


(2.5) 


O 

,  Vrx* 

Rx  *  Vry*Rv  *  Vrz*Rz 

(2.6) 

J\ 

~  '  y  ^ 

(*x  +  Ry  + 

Where  V  * 

r  a 

Vtx  ‘ 

Vmx 

V 

ry 

« 

X 

4-> 

> 

V 

my 

V  ,  * 
rz 

Vtz  * 

V 

mz 

are  the  relative  velocity  components  in  the  inertial  frame. 

In  the  geometry,  the  relative  angles  ^  and  also 

will  be  expressed  by  the  relative  vectors: 


C-7) 


,  ,  -1  Rv 

*r  -  -tan  r 


Where  0R  and  are  the  elevation  and  the  azimuth  angles 

of  the  relative  sight  line  vector  in  the  inertial  frame. 

Then  the  rate  of  change  of  both  angles  are 


i  .?n 

*R  R 


-  sinevcos-4/vVrx  -  sin6vsin<|,vVry  +  c0seyVrz 

(RXZ  * 


(2.8) 


r 'p  _ 


*sin:V'rx  +  cos*vVi 


RcosS.  (R  *  +  R  R  " 

V  v  X  y  +  Z  )C0S< 


(2-9) 


Where  VRQ  and  VR  are  the  rates  of  change  in  ey  and 
ipv  components,  is  shown  in  Figure  2.3.  These  quantities 
will  be  used  later  in  developing  the  estimator  state 

equations  and  will  also  be  used  to  represent  the  variables 

* 

for  the  missile  seeker.  The  quantities  ®r  >  6r  >  ^r  »  ’J'r  » 
R  and  R  (corrupted  by  noise)  are  the  set  of  measurements 
assumed  available  from  the  missile  seeker. 


D.  KINEMATICS  OF  CONTROL  VECTOR 

Figures  2.4  and  2.5  are  views  from  near  and  side  of 


missile.  a^  and  $  represent  the  magnitudes  of  the 


21 


DEROLLED 

-Z-AXIS 


DEROLLED  Y-AXIS 


FIGURE  2.4  KINEMATICS  OF  CONTROL  VECTORS  VIEWED 
FROM  REAR  OF  MISSILE. 


FIGURE  2.5  KINEMATICS  OF  CONTROL  VECTORS  VIEWED 
FROM  SIDE  OF  MISSILE. 


missile  acceleration  due  to  the  acceleration  command  (ac) 
and  the  bank  angle  due  to  the  roll  rate  command  (Pc)  each 
other,  g  represents  the  gravitational  effect  at  the  center 
of  gravity  of  the  missile.  As  mentioned  in  control  method, 
the  sign  of  a  is  such  that  it  is  positive  upward  with 
respect  to  the  wings.  Figure  2.5  shows  that  the  accelera¬ 
tion  normal  to  the  velocity  vector  has  the  components  due 
to  the  elevator  control  input  and  the  effect  of  the  gravity. 
The  total  magnitude  of  the  acceleration  normal  to  the 
velocity  vector  can  be  represented  as 

a  *cos4>  -  g*sin6 
m  45  v 

The  velocity  and  acceleration  components  for  simple 
circular  motion  can  be  represented  as 

an  *  v2/P 

See  Reference  1. 

Using  this  definition,  the  differential  equations 
describing  the  rate  of  change  of  the  flight  path  angles  are 


■ 


k 


am  SiniJ) 

J,  =  JB _ 

Vv  V  cose 
m  v 


C2.ll) 


iNote  am*cos<{)  -  g*cos9v  is  the  magnitude  of  the  normal 
acceleration  component  to  the  velocity  vector  in  missile 
side  plane  and  am*sin<f  is  the  magnitude  of  the  normal 
acceleration  component  to  the  velocity  vector  in  the  x-y 
plane.  The  acceleration  component  of  the  missile  velocity 
V  has  only  the  component  of  the  effect  of  gravity  in  the 
missile  side  plane  under  the  assumption  that  the  missile 
velocity  is  constant.  The  differential  equation  is: 

Vm  =  -g*sin9v  (2.12) 


The  rates  of  change  of  the  missile  position  components 
(Xm,  Ym,  2m)  can  be  obtained  by  integrating  the  components 
of  missile  velocity.  From  Figure  2.2,  the  components  of 


missile  velocity  (Vmx, 

Vmy,  Ymz)  in  inertial  frame 

are 

given  by 

Xm  = 

V 

=  V  *COS0  *COS’)J 

(2.13) 

m 

rax 

m  v  rv 

Y  = 

V 

=  V  *cos0  *sin^ 

(2.14) 

m 

my 

m  v  rv 

2  = 

Y 

=  -  V  *sin9 

(2.15) 

m 

mz 

m  v 

Let  the  response 

of 

the  missile  to  input  command 

in 

normal  acceleration  and  roll  rate  be  considered.  If  a 
demand  is  made  on  a  missile  for  a  transverse  acceleration, 
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it  is  initiated  by  sending  a  signal  to  the  appropriate  con¬ 


trol  surface.  The  missile  acceleration  a  follows  the 

m 

demanded  acceleration  a^  in  a  manner  which  may  be 
characterized  by  a  frequency  (weather  cock  frequency) 

and  a  damping  factor  P  .  For  simplicity,  in  present 
work,  only  the  first  order  lags  with  time  constant  and 

tq  are  assumed,  the  equation  of  the  first-order  system  can 
be  represented  as  follows: 

X(t)  +  l./x*X(t))  =  f(t) 

Taking  the  Laplace  transformation  on  has 

X(S)  =  F(S)  /  (S  +  l .  /T) 


with  zero  initial  condition.  Applying  these  to  the  BTT 
missile  control  system. 


T 


+ 


m 


*  a. 


(2.16) 


T  *P  +  P  *  P 
P  C 


(2.17) 


Rearranging  the  differential  equation. 


Where  a„  and  P 
c  c 


are  the  acceleration  and  roll  rate 


commands.  Equation  (2.20)  is  by  definition.  From  the 
above  differential  equation?  the  control  ratios  are 
transformed  as  follows: 


%,(=:> _ i_ 

a  (S)  t  S  +  1 
c v  a 


(2.21) 


P  (S)  _  1 

pTTsT  ~  ips  +  i 


(2.22) 


Equations  (2.10)  through  (2.22)  constitutes  the  dynamic 
equations  describing  missile  motion  in  response  to  the 
command  inputs  ac  and  Pc  .  The  solution  of  these 
equations  provide  missile  position  (Xm  ,  Ym  ,  Zm)  , 
orientation  and  magnitude  of  velocity  ,  0y  ,  Vm)  and 

orientation  and  magnitude  of  control  acceleration  ($  ,  am) 


E.  CO-ORDINATE  TRANSFORMATION 

Given  that  the  missile  is  on  guided  flight  to  compute 
the  magnitudes  of  the  control  inputs  at  any  moment,  it  is 
necessary  to  transform  the  instantaneous  values  from  the 
sensors  in  the  inertial  frame  to  the  instantaneous  values  in 
body  frame.  The  transformation  from  one  frame  to  another 
can  be  accomplished  through  a  transformation  matrix.  The 
transformation  matrix  wiill  be  developed. 

Define  an  inertial  coordinate  system  with  unit  vectors 
I,  J,  K.  Also, define  a  missile  body  axis  coordinate  system 
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with  unit  vectors  i,  j,  k.  It  is  desired  to  find  the  trans¬ 
formation  between  the  I,  J,  K  system  and  the  i,  j,  k  system. 

The  I,  J,  K  system  can  be  thought  to  be  oriented  so 
that  I  points  north,  J  east,  and  K  down.  Similarly,  the 
i,  j,  k  system  has  i  along  the  longitudinal  axis  of  the 
missile  (which  by  assumption  is  along  the  velocity  vector 
Vm) ,  j  out  the  right  wing,  and  k  down. 

Consider  three  successive  rotations  ,  9y  and  i>  . 

The  tj/  rotation  is  about  the  inertial  z  axis,  and 
transforms  to  an  intermediate  axis  systems  i^,  j^,  k^. 

The  0V  rotation  is  about  the  axis  and  transforms  to  an 
axis  system  i7,  j^,  k7.  The  last  rotation  <p  about  the 
i7  axis  transforms  from  flight  path  axes  to  missile  body 
axes.  As  a  result,  the  total  transformation  from  inertial 
to  body  axes  is  given  by 


[$J[ev][^v3  j  -  [A]  j 


(2.23) 


Where  A  is  the  total  transformation  matrix  from  inertial  to 
body  frame.  The  elements  of  the  a  matrix  from  Reference  1 


all  *  cos (Q  )  *cos  (i|>  ) 
al2  =  cos  ( 0  )  *sin(ifi  ) 
al3  *  -sin(6v) 


III.  AN  OPTIMAL  CONTROLLER 


Although  this  work  is  primarily  concerned  with  the 
effect  of  an  optimal  estimator  on  the  performance  of  a  bank 
to-turn  missile,  we  will  first  outline  the  theory  of  an 
optimal  controller  in  this  section.  The  guidance  laws 
considered  here  were  first  developed  by  Stallard  [Ref.  10]. 
The  following  is  a  brief  outline  of  his  work. 

A.  GEOMETRY 

Figure  3.1  shows  the  applicable  geometry  representing 
the  Proj ected-Zero-Control  miss  (PZC  miss)  distance  in  the 
terminal  state.  The  problem  is  considered  to  be  three- 
dimensional  and  two  major  axes  system  are  used: 

1.  A  seeker-oriented  axis  system  for  estimation  or 
measurement,  and 

2.  The  three  principal  axes  of  the  missile  for  the 
control  problem, 

Xb ,  Yb,  Zb  denote  the  target  coordinates  relative  to  the 
missile  along  its  principal  axes.  The  body  system  is  repre 
sented  a  time  T  *  0  by  a  set  of  Eulian  angles  which 
relate  it  to  an  earth  fixed  system.  The  angle  Af  is  the 
incremental  roll  angle  from  the  present  missile  axes  to  any 
future  orientation  at  time  T  ,  shown  dashed  in  Figure  3.1. 
The  PZC  miss  represents  the  total  Pro j ected- Zero  effort  - 
miss  distance  and  is  defined  as  the  miss  distance  which 


would  occur  if  no  further  control  was  exercised  on  the 
missile  at  T  =  Ti  .  It  is  composed  of  a  Y-components 
(P2CY  miss)  and  a  2-component  (P2C2  miss).  Under  the 
assumption  of  constant  velocity  the  X-component  (PZCX  miss) 
is  represented  by  the  distance  that  the  missile  has  to  fly 
at  any  time,  since  it  is  not  effected  by  control,  it  is 
not  shown  in  Figure  3.1. 


n 


FIGURE  3.1  COMPOXEXTS  OF  P2C ,  MISS  DISTANCE. 
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B.  SUMMARIZED  DESCRIPTION’  OF  AN  OPTIMAL  CONTROLLER 


Given  a  tine-varying  linear  systen  of  the  form, 

X  =  F (t) X  +  G(t)U  (3.1) 


Where  X  is  n-conponent  state  vector 

U  is  m-coraponent  control  vector. 

An  optimal  control  is  one  that  minimizes  a  performance 
index  made  up  of  a  quadratic  form  in  the  terminal  state 
plus  an  integral  form  of  quadratic  forms  in  the  control  and 
the  state: 


J  = 


I  (xTs 


f* 


T  = 


Ti) 


UTBU)  dt 


(3.2) 


Where  the  S^.  and  A  are  positive  semidefinite  matrix  and 
B  is  a  positive  definite  matrix.  An  appropriate  choice  of 
these  matrices  must  be  made  to  obtain  acceptable  level  of 
X(Ti)  ,  (X(t)  and  U(t)  . 

The  guidance  law  determined  in  Reference  8  is  of  the 
fo  rm : 


uc  ■  (ac  -  pDT 

With  the  performance  index  chosen  (the  weighting  on  the 
state  vector  was  not  considered  in  Reference  8. 


Where  the  first  term  is  one-half  the  required  miss  distance, 
the  second  term  is  included  in  order  to  eliminate  infinite 


energy  solutions.  B  is  a  2  by  2  positive  definite, 
symmetric  weighting  matrix  and  of  the  form 


B 


B1  0 
0  B2 


(3.4) 


with  1/Bjj  =  (Ti  -  To)X  maximum  acceptable  value  of 
2 

|Ui(T) !  .  A  detailed  derivative  of  the  optimal  guidance 

law  is  given  in  Reference  10. 

The  solutions  of  the  above  equation  from  Reference  8 
give  the  optimal  acceleration  command  at  time  T  between 
TO  and  Ti  as 


3T. 
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3B,  +  T 

1  go 


I  C-  Mbz> 


(3.5) 


Where  Tg0  is  the  remaining  time  to  go  until  intercept  and 
Mfa2  is  the  PZC-miss  along  the  missile  normal  force  axis  and 
the  negative  sign  is  due  to  definition  of  the  normal 


The  optimal  controller  is  implemented  in  the  present  work 
using  Equations  (3.5),  (5.7),  and  (3.9).  In  computing  the 
control  input,  Tao  ,  Y  and  Z-component  of  PIC-miss  and 
the  elements  of  the  weighting  matrix  are  needed. 

C.  APPROXIMATING  TIME-TO-GO 

It  is  necessary  for  the  Time-To-Go  and  the  Projected- 
Miss - Dis tance  in  missile  body  axes  to  be  specified  at  time 
T  to  compute  the  optimum  control  input.  The  Time-To-Go 
will  be  defined  at  first,  then  the  state  variables  will  be 
introduced  as  means  for  formulating  the  necessary  equations. 
The  target  acceleration  components  of  the  state  variables 
will  be  modeled  since  it  is  assumed  that  any  information 
about  the  target  acceleration  cannot  be  obtained  from  the 
active  seeker  of  the  missile.  Then  using  this  target 
acceleration  model  the  algorithm  for  computing  the  other 
state  variables  will  be  developed. 

The  Time-To-Go  until  intercept,  T^  ,  is  required  for 
computing  the  control  input.  Let  the  time  to  intercept  be 
defined  as  the  time  to  minimum  range.  Under  the  assumption 
that  the  relative  velocity  vectors  may  be  considered 
constant  at  its  present  value  until  intercept.  T  can  be 
found  at  time,  T  ,  as: 


"V  J 


Where  R  is  the  present  position  of  the  target  relative  to 
the  missile  and  the  factor  in  brackets  in  the  component  of 
V  along  R  .  An  alternate  form  of  T  is: 


Ti  = 


Vrx  +  V'ry  +  V'rz 

Y  +  V  "  +  V  c\ 
rx  ry  rz  1 


(5.11) 


The  negative  sign  comes  from  the  sign  of  V  as  shown  in 
Figure  2.2  and  Equation  (2.6). 

D.  PROJECTED- ZERO- CONTROL  MISS  ( P2C  MISS)  DISTANCE 

The  Proj ected- Zero-Control  miss  distance  (PZC-miss)  is 
defined  as  the  miss  distance  which  would  occur  if  there 
were  no  further  control  input  after  time  T  . 

In  this  work  the  state  vectors  are  defined  as  the 
following  state  variables 

X  =  (Rx,  Ry,  Rz,  VrJ£,  Vry,  Vr_,  atx,  aty ,  atz)T 

Where  R  is  the  relative  missile  distance  to  target 

V  is  the  relative  missile  velocity  to  target 

is  the  target  acceleration 

A  model  for  target  acceleration  has  been  assumed.  This 
is  given  as  follows: 


at(T)  =  atoe’  a 


a  >  0 


(3.12) 
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Where  a  is  tne  intial  target  acceleration  and  a  is  the 
reciprocal  of  the  acceleration  time  constant.  The  initial 
target  acceleration,  4g's,  will  be  assumed  for  each 
scenario,  and  a  =  1/20  for  evasive  maneuver  as  indicated 
in  Reference  10  will  be  chosen.  These  values  will  be  used 
in  the  simulation  of  the  optimal  controller. 

With  the  assumption  that  a  state  estimator  is  available 
to  provide  estimates  of  the  state  vector  at  current  time  (t) 
the  values  of  the  state  variables  for  no  control  inputs  at 
future  time,  can  be  obtained  easily  using  the  integration, 


1 .  e .  , 


V_(t)  =  V  (tO)  +J  a  (f) dx 
r  r  to  Z 


R(t)  «  R(t0)  ♦/  V  (t  )  dt 

+■  n  * 


The  results  of  these  integrations  are 


lr(t)  -  Vr(t0)  ♦  1/a  (1  -  e  ■a[t  ’  t0))at(t0)  (3.13) 


R(t)  =  R(t0)  +  (t  -  tO)*V  (to) 


(3.14) 


[  1  /  a 2  ( :t  ( t  -  tO)  -  1  +  e  'a(t  '  t0))|*at(t( 


Where  tO  and  t  denote  current  time  and  future  time.  It 
should  be  noticed  that  above  integration  is  based  on  the 
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assumption  of  a  single  physical  dimension.  If  this  assump¬ 
tion  is  to  be  applied  to  three-dimensional  problem,  it  must 
be  concluded  that  there  are  no  coupling  among  the  X,  Y,  Z 
dimension.  Let  T  =  t  -  t  ,  then  the  components  of  the 
PZC-miss  in  inertial  axes  at  current  time  (tO)  can  be 
expressed  as 
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(3.15) 


0 

0 

(l/2)*g*T 


go  J 


Where  f  =  f  =  f  =  l/a^[aT  -  1  +  e  a^go]  and  the  last 

X  >  go 

term  accounts  for  the  free  fall  effect  of  gravity  in  the 
vertical  axis . 

The  components  of  the  PZC-miss  in  missile  body  axes  for 


control  input  calculation  are  obtained  via  the  Euler  trans¬ 
formation  matrix: 
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Where  it  is  assumed  that  the  missile  has  exact  knowledge  of 
its  orientation  angles  4>  ,  3  ,  $  ,  from  the  accurate 

V  v 

inertial  measurement  units. 

The  values  of  optimal  control  vectors  (a  and  P  )  , 
r  c  c 

then,  are  calculated  from  time  to  go  (T  )  PZC-miss  in 

°  OO 

o 

body  axes  and  weighting  matrix  components  (B1  and  B2)  which 
may  be  determined  according  to  the  desired  miss  distance  and 
intercepting  time.  In  a  later  section,  the  weighting 
matrix  will  be  developed  and  the  scenarios  will  be  set  up 
for  computer  simulation. 

E.  SIMULATION  RESULTS 

The  guidance  law  in  this  work  was  tested  on  a  simple 
six- degree  -  of - freedom  simulation  using  only  plausible 
numbers  without  reference  to  any  actual  missile  or  design. 

Before  the  computer  simulations,  the  weighting  matrix 
of  the  performance  index  selected  in  Reference  8  will  be 
calculated  using  the  initial  state,  .\(0)  ,  the  missile 
operating  time  imposed,  T  (5  sec)  and  the  constraint  input 
vectors.  It  is  assumed  that  the  missile  has  a  velocity 
advantage  of  two  over  that  of  the  target,  a  zero-  and  0.5 
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second-time  lag  auto  pilots  for  pitch  with  limits  of  -3g's, 
+20g's  and  for  roll  rate  with  limits  of  +6.28  rad/sec. 

The  allowable  miss  distance  which  is  usually  determined 
by  sizing  of  the  warhead  was  chosen  as  5  meters.  Which 
makes  the  first  term  of  the  performance  index  Equation 
(3.3)  equal  to  12.5  meters.  The  acceleration  term  and  the 
roll-rate  term  in  the  performance  index  are  set  equal  to 
this  respectively  (one-half  of  square  mean  miss  distance) , 
then 

(a  -  a,)2  dt  =  12.5  m2  (3.16) 

v  c  d 

0 

(Pc2)  dt  =  12.5  m2  (3.17) 

0 

The  nominal  time  of  flight  for  the  terminal  phase  is  set 

equal  to  5  seconds.  The  components  of  the  weighting 

matrix  B  were  based  on  the  following  assumptions.  It  is 

2 

assumed  for  (4g's)  as  a  mean-square  value  of  (ac  -  a^)  and 
2 

(2  rad/sec)  as  a  mean-square  roll  rate  to  be  acceptable. 
Then  the  components  of  weighting  matrix  are: 

B1  =  0.003254  sec3 
2 

B  *  1.25m  sec 
2 
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The  commands  a  and  P  are  computed  every  0.05  sec 

and  are  frozen  when  the  Time-To-Go  falls  below  this  interval. 

In  order  to  simulate  the  optimal  controller  appropriate 
scenarios  were  developed  below  with  restriction  that  the 
Time-To-Go  was  5  sec  and  the  speed  of  missile  is  two  times 
faster  than  that  of  target.  Three  scenarios;  tail-chase, 
head-on,  s ide- approach  engagements,  were  chosen  as  possible 
cases  in  this  work. 

In  case  1,  tail-chase  engagement,  the  target  flies  100m 
above  the  missile  and  2500m  uprange  with  a  horizontal  east¬ 
ward  acceleration  of  4g's.  Figure  5.3  shows  the  case  2. 

In  case  5,  s ide- approach  engagement,  the  target  flies  400m 
above  the  missile  and  1900m  uprange  with  a  horizontal  north¬ 
ward  acceleration  of  4g's  at  t  =  0  as  shown  in  Figure  3.4. 
For  simplicity,  when  the  scenarios  were  set  up,  lOOOm/s  was 
taken  as  a  missile  speed,  500m/s  as  a  target  speed  in  case 
1  and  case  2,  600m/s  as  a  missile  speed,  300m/s  as  a  target 
speed  in  case  3,  and  4g's  as  a  target  acceleration  for  all 
three  cases.  When  simulation  was  executing,  this  first 
simulation  assumed  perfect  state  feedback  without  the  use  of 
an  optimal  estimator. 

The  simulation  results  are  shown  in  Figure  3.5  through 
Figure  3.10.  At  first,  the  optimal  controller  with  no  time 
lag  was  tested  for  three  scenarios,  as  shown  in  Figures  3.5 
and  3.6.  For  scenario  1,  the  normal  acceleration  command 
in  initial  phase  is  very  large,  then  decreases  to  zero  at 
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FIGURE  3-5.  ROLL  RATE  COMMAND  VS  TIME 


at  the  final  time.  For  scenario  2,  the  acceleration  is  zero 
in  the  initial  phase,  then  increases  and  converges  to  zero 
finally.  For  scenario  5,  the  acceleration  starts  from  zero, 
but,  increases  rapidly  up  to  the  maximum  acceleration  then 
freezes  for  seconds  before  decreasing.  The  acceleration 
does  not  converges  to  zero  at  the  final  time.  Comparing 
these  results  with  the  variation  of  the  roll  rate  command 
in  Figure  3.6,  the  opposite  trends  are  observed  in  the 
results  of  scenarios  1  and  2,  i.e.,  if  an  acceleration 
command  in  the  initial  phase  is  large,  a  roll  rate  command 
is  small.  The  acceleration  command  in  Equation  (3.5)  is  a 
linear  function  of  PZC-miss  in  of  z-direction  of  the  missile 
coordinate  system.  For  scenario  1,  PAC-miss  in  the 
z-direction  is  relatively  large,  this  causes  a  large  accel¬ 
eration  command  and  a  small  roll  rate  command  in  initial 
phase.  For  scenario  2,  the  relative  magnitude  of  the  PZC- 
miss  in  the  z-direction  is  small.  This  cause  a  large  roll 
rate  command  in  the  intial  phase.  For  scenario  5,  the 
PZC-miss  in  v-direction  is  of  the  same  magnitude  as  the 
relative  distance  as  the  Time-To-Go  decreases.  This  means 
that  a  larger  acceleration  command  than  the  maximum  value 
allowable  in  the  missile  body  coordinate  system  is  required 
at  some  of  instance  time  during  the  missile  flying.  Thus 
the  acceleration  does  not  converge  to  zero  at  the  final  time. 
In  the  derivation  of  the  equations  for  control  commands,  a 
small  angle  approximation  was  assumed.  However,  in  the 


scenario  3  this  assumption  is  questionable.  The  roll  rate 
command  computed  as  shown  in  Figure  3.6  shows  a  sine  wave 
characteristics  which  is  possibly  due  to  overcorrection. 

The  simulation  results  for  a  0.5  second  time  lage  are 
shown  in  Figure  5.8  and  Figure  3.9.  The  time  lag  causes 
much  larger  acceleration  control  commands.  This  includes 
overcorrections,  so  that  the  control  commands  do  not 
converge  in  the  final  phase.  The  other  characteristics 
follow  the  explanations  given  in  the  analysis  of  the  cases 
with  no  time  lag. 

For  all  three  scenarios  the  miss  distances  were  obtained 
in  the  simulation  of  both  cases,  no-  and  0.5  second-time 
lag.  The  miss  distances  are  below  0.5m  for  all  three 
scenarios  in  the  case  of  no-lag.  In  the  no  lag  case  the 
control  laws  gave  excellent  results.  In  the  lag  case,  the 
miss  distances  are  5.9m  for  tail-chase  engagement,  0.8m  for 
head-on  engagement,  and  32.9m  for  side-approach  engagement. 
The  miss  distance  of  the  side- approach  engagement  with  a 
different  maximum  acceleration  limit  fo  23g's  results  in  a 
miss  distance  of  less  than  5m.  This  higher  maximum  acceler¬ 
ation  limit  will  be  used  in  the  scenario  3  during  the  Simula 
tions  for  testing  the  performances  of  whole  control  system 
in  Chapter  5. 

The  previous  results  are  summarised  as  follows:  (1)  the 
tracking  of  a  target  is  very  dependent  of  the  missile  and 
target  geometry,  (2)  the  capabilities  of  the  missile 
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(maximum  acceleartion  and  maximum  roll  rate  limits  and  time 
constant)  in  maneuvering  against  the  target  is  very  important 
factors,  (5)  higher  maximum  accelerations  are  required  in 
any  side  approach  engagement,  and  (4)  the  assumption  of 
small  in  the  optimal  control  law  is  a  limitation  in  the 

general  case.  Nevertheless,  the  control  laws  developed 
from  the  optimal  control  theory  for  both  the  no  lag  and  the 
0.5  second-time  lag  cases  yielded  successful  results  for 
the  specified  scenarios.  Since  the  optimal  controller  with 
0.5  second-time  lage  auto  pilots  resulted  in  a  miss  distance 
of  approximately  5m  with  complete  state  variables  feedback, 
it  is  expected  that  the  inclusion  of  an  optimal  estimator  for 
these  particular  scenarios  will  results  in  still  larger  miss 
dis  tances . 
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IV.  OPTIMAL  ESTIMATOR 


A.  DESCRIPTION 

Without  full  state  variable  feedback  it  is  necessary  to 
estimate  the  state  variables.  This  chapter  describes  the 
Kalman  filters  used  in  the  state  estimation.  These  filters 
were  designed  for  respectively:  (a)  two  measured  states, 
(b)  four  measured  states,  and  (c)  six  measured  states. 

All  measurements  have  noise  superimposed  on  them.  It  i 
assumed  that  the  nature  of  the  noise  source  is  known 
completely . 

The  Kalman  filter  is  the  optimal  estimator  of  the 
states.  Since  the  measurements  are  non-linear,  it  was 
necessary  to  use  the  extended  Kalman  filter  theory.  A 
first  order  extended  filter  was  assumed  in  this  work.  In 
developing  the  filter  algorithm,  the  non-linear  measurement 
equation  is  linearized  about  the  most  recent  state  estimate 
and  then  the  Kalman  filter  algorithm  is  applied.  This  is  a 
extended  Kalman  filter  algorithm.  The  results  of  the 
extended  Kalman  filter  will  yield  suboptimal  state  estima¬ 
tion. 

In  practice,  the  missile  will  have  sensors  of  various 
types  that  provide  inputs  which  are  related  in  some  fashion 
to  the  states  to  be  estimated.  The  main  emphasis  in  the 
present  work  will  be  on  active  seeker  which  is  assumed  to 


provide  the  optimal  estimator  with  measurements  of  line-of- 


sight  angles  for  the  case  of  the  missile  having  two  measure¬ 
ment  sensors.  For  the  case  of  four  measurements  one  lias  two 
of  line-of-sight  angles,  range  and  time  rate  change  of 
range.  For  the  case  of  six-measurements  one  has  two  of 
line-of-sight  angles,  two  of  time  rate  change  of  the  line- 
of-sight  angles,  range  and  its  time  rate  change.  When  the 
measurement  equations  and  the  Jacobean  matrix  are  developed 
the  six -measurement  case  will  be  derived  since  the  other 
cases  have  the  same  kind  but  less  number  of  measurement 
sensors . 

We  will  first  consider  the  state  equations  and  then  the 
measurement  equations.  This  will  be  followed  by  a  discus¬ 
sion  of  the  extended  Kalman  filter.  Measurement  error 
effects  on  the  initialization  of  the  Kalman  filter  will  be 
analyzed.  Finally,  the  estimator  for  the  measurement  cases 
will  be  simulated  in  order  to  test  the  effects  of  the  imple¬ 
mentation  of  the  different  measurment  vectors  on  the 
extended  Kalman  filter. 


B.  STATE  EQUATIONS 

The  state  vector  is  limited  for  convenience  to  the 
following  state  variables: 


(Rx  ' 


Rz  ’  ^rx’  ^ry’  ^  rz  ’  atx’  aty’  atZ" 


(4.1) 
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Where  ,  X?  =  R^  ,  X-  =  R,  are  the  components  of 

relative  position,  X.  =  V  ,  Xr  =  V  ,  X,  =  V  the  com- 
1  *  4  rx  5  ry  6  rz 

ponents  of  relative  velocity,  and  X_  =  a.  ,  X0  =  a.  , 

/  lA  O  L  } 

Xg  =  a^,  the  components  of  target  acceleration. 

In  order  to  develop  the  state  equations,  dynamic  equa¬ 
tions  of  target  motion  are  needed.  The  target  model 
selected  for  tracking  applications  must  be  sufficiently 
simple  to  permit  ready  implementation  in  weapons  system  fo 
which  computation  time  is  at  a  premium  yet  sufficiently 
sophisticated  to  provide  satisfactory  tracking  accuracy. 

To  meet  these  requirements,  the  model  described  by  Singer 
[Ref.  10]  will  be  used  in  this  work.  The  model  will  be 
presented  for  a  single  spatial,  dimension  in  order  to 
enable  accurate  tracking  performance  estimates  to  be  made 
for  a  variety  of  sensor  measurement.  If  the  targets  under 
consideration  normally  move  at  constant  velocity,  the 
accelerations  due  to  turns  or  evasive  maneuvers  may  be 
viewed  as  perturbations  upon  the  constant  velocity  trajec¬ 
tory.  The  target  acceleration  a(t)  therefore  will  be 
termed  the  target  maneuver  variable,  in  a  single  physical 
dimension.  The  target  maneuver  capability  can  be  satis¬ 
factorily  specified  by  two  quantities:  the  variance  or 
magnitude  of  the  target  maneuver  and  the  time  constant, 
or  duration,  of  the  target  maneuver.  Hence  the  target 
acceleration,  namely,  the  target  maneuver  is  correlated 


in  time.  A  typical  representative  model  of  the  target 
acceleration  as  a  correlation  function,  can  be  expressed 
as  follows : 

at(x)  =  E[at(t)at(t  +  x)]  =  om"e  a  >  0  (4.2) 

Where  the  a  ^  is  the  variance  of  the  target  acceleration 
m 

and  a  is  the  reciprocal  of  the  maneuver  time  constant. 

The  variance  o  “  of  the  model  will  be  approximated  using 
m 

the  following  formula  which  represents  the  acceleration 
probability  density  suggested  by  Singer  [Ref.  10]: 

°m2  '  [l  *  4P»ax  •  Pol  '4-3> 

where  a  is  a  maximum  acceleration  rate,  Pmax  is  the 
insx 

probability  of  maneuvering  at  amax  ,  and  PQ  is  the 
probability  of  no  target  maneuvering.  And  as  mention  in 
previous  chapter,  a  -  1/20  for  an  evasive  maneuver  will 
be  used  in  this  work. 

The  above  process  can  be  modeled  as  the  output  of  a  low- 
pass  filter  driven  by  white  noise  and  it  can  be  described 


by  the  equation 


.  *.  *W  *- 


i 


where  c  .  (t)  ,  the  correlation  function  of  the  white  noise 
w 

input  satisfies 


aw“(T)  =  2acm‘'5(r) 


(4.5) 


The  acceleration  model  above  will  be  applied  to  each 
X,  Y,  Z  dimension  with  no  cross  coupling,  the  system  state 
equations  can  be  written  then  as: 

X-FX+a  +  W 


(4.7) 

The  Equation  (4.7)  is  a  linear  state  equation  with  the 
assumption  that  the  forcing  function  a  vector,  i.e.,  missile 
acceleration,  can  be  precisely  measured  and  also  resolved 
into  X,  Y,  Z  components  in  the  inertial  frame. 
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Many  sensors  have  a  constant  data  rate,  sampling  the 
data  every  T  seconds.  If  the  above  system  equation  is 
represented  in  discrete  equations.  The  appropriate  system 
equations  of  motion  may  be  expressed  by 


X(K  +  1)  =  $ (T) X(K)  +  B(K)  +  U(k') 


(4.8) 


where  $(T)  is  the  target  state  transition  matrix,  3(K) 
the  deterministic  forcing  input  vector  which  is  composed  of 
first  and  second  integration  and  U(K)  the  inhomogeneous 
driving  input  due  to  white  noise.  Since  it  is  assumed  that 
the  missile  acceleration  is  a  known  deterministic  forcing 
function,  the  system  equation  of  motion  can  be  obtained  by 
direct  integration  in  each  single  dimension  with  E|K(t)|=0 
The  desired  integration  in  discrete  form  yields 


R(tO+T) 


1  T  l/a2(-l  +  aT  +  e’aT)  I  |R(tO) 


|Vr(tO+T)  =  0  1 


l/a(l-e’aT) 


vr(to) 


aT(tO+T) 


at(tO) 


-xr°„:T 


tOam  (t)dtdt 


f  tO+T 

J  a  (t)dt 


X(K  +  1)  =  (t  ,  a  )  X(K)  +  am(K)  (4-9) 

Where  K  +  1  =  (K  +  1)T  =  t  +  T,  K  =  KT  =  tO  . 


.s'ext,  if  the  driving  input  U(t)  vector  is  considered, 
this  input  is  not  a  sampled  version  of  the  continuous  time 
white  noise  input  W(t)  .  Since  W(t)  is  white  noise, 

U(k)  also  should  be  a  discrete  time  white  noise  sequence, 
that  is,  E [U(K) U(K+i) ]  =  0  for  i  ^  0  .  Then  as  shown 
in  Reference  10,  U  and  W  are  related  by 

(K+l) 

U(k)  -J  $| [K  +  l)t  -  t |  W(t)dt  (4.10) 

KT 

Where  U(K)  is  a  white  noise  sequence  with  covariance 
mat  rix 

E[U(K)U(K)T]  =  Q  (4.11) 


Following  the  above  derivation,  it  can  be  easily 
expanded  to  the  case  of  three  independent  dimensions  in  X  , 
Y,  Z.  For  otx  =  ay  =  a,  =  a  ,  and  ox  =  ay  =  o z  =  a 
since  the  model  of  target  acceleration  is  for  one  dimension, 
the  transition  matrix  is 


I  TI  f2I 


1  f2X 


(4.12) 
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Where  I  and  0  are  3  bv  5  identify  and  null  mat r 


f:  =  l/a~ (aT  -  1  +  e’aT) 


f,  -  1/ot  ( 1  -  e'aT) 


f  3  =  e 


-aT 


The  forcing  acceleration  vector  is 
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The  covariance  matrix  of  white  noise  U(K)  is 


Whe  re 
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C.  MEASUREMENT  EQUATION 

In  the  case  of  the  estimator  with  six -measurement  sen¬ 
sors  the  present  study  assumes  that  the  active  seekers  of 
missile  give  the  measurements  of  the  sightline  angles 


,  the  time  rate  change  of  the  sight  line  angles 
9^  ,  ii'  ,  the  relative  range  R  ,  and  the  time  rate  change 
of  the  relative  range  R  contaminated  'ey  the  Gaussian  white 
noise.  Then  the  equations  of  measurements  in  discrete  time 
form  are 

2 (K)  =  h(x,K)  +  V(K)  (4.15) 

where  Z_(K)  is  a  sequence  of  measurement  vector  contaminated 
by  white  noise  V(K)  and 

elevation  angle  6^ 

elevation  angle  rate  0r 

azimuth  angle  y 

h(x)  =  -  -r 

azimuth  angle  rate  4>r 

range  R 

range  rate  R 

The  missile  seekers  actually  measure  the  data  about  the 
seeker  axes.  Thus,  throughout  this  development  the  assump¬ 
tion  is  made  that  the  missile  possesses  an  inertial  measure¬ 
ment  unit  that  accurately  specifies  the  missile's  orienta¬ 
tion  in  space  so  that  a  transformation  from  seeker  to 
inertial  coordinates  can  be  made. 

Referring  to  Figures  2.2  and  2.3,  the  six-measurement 
vectors  can  be  written 


I 


The  measurement  vector  can  be  expressed  in  terms  of  stati 
vector 


j-  tan 
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(4.16) 


The  measurement  noise  covariance  matrix  for  the  active 
seeker  can  be  written 


Og  0  0  0  0  0 

0  Og2  0  0  0  0 

0  0  a ,  2  0  0  0 

0  0  0  s’.2  0  0 

0  0  0  0  a  ^2  0 

0  0  0  0  0  a , 


E(V(K)V(K)Tj  =  R  = 


(4.17) 


1  *  "  ' 


Where  the  diagonal  elements  are  the  variances  of  th.e  indi¬ 
vidual  measured  quantities. 

As  shown  in  Equation  (4.16),  the  measurement  vectors  are 
the  non-linear  functions  of  state  vectors.  It  is  impractical 
to  implement  in  non-linear  form  because  the  computation  of 
gain  K  and  error  covariance  matrix  P  in  update  process 
in  extended  Kalman  filter  algorithm  is  not  possible.  To 
simplify  this  computation,  and  to  implement  the  extended 
Kalman  filter,  the  Jacobean  or  matrix  of  partial  derivatives 
H  will  be  determined,  i.e., 


(4.18) 


The  Jacobean  will  be  a  6  by  9  matrix  because  h 
with  six  elements  and  A  is  a  vector  with  nine 
Performing  the  operation,  the  components  of  the 
are 
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D.  THE  EXTENDED  KALMAN  FILTER 


Optimal  estimators  that  minimise  the  estimation  error, 
can  be  divided  into  two  processes,  the  one  is  the  updating 
process  to  estimate  the  state  vector  at  the  current  time, 
based  upon  all  past  measurements,  the  other  is  the  extra¬ 
polating  process  to  estimate  the  state  at  some  future  time. 
Optimal  estimation  procedure  for  linear  state  equations  and 
non-linear  measurement  in  the  form  of  discrete  time  will  be 
described  referring  to  Reference  6.  Figure  4.1,  a  timing 
diagram  shows  the  flow  of  the  various  quantities  involved 
in  the  discrete  optimal  filter  equations. 

The  discrete  system  equation  whose  state  at  time  kt  is 
denoted  by  simply  X(K)  ,  where  U(K)  is  a  zero  mean, 
white  sequence  of  covariance  Q(K)  ,  is 

X(K)  -  i  X(K  -  1) ,  +  B(K  -  1)  ♦  U(K  -  1)  (4.21) 

where  4>  is  a  transition  matrix 

E[U(K)]  =  0  (4.22) 

E[U(i)  U(  j)T]  =  Q5  i  j  (4.23) 

The  measurements  are  the  non-linear  function  of  the  system 
state  variables,  corrupted  by  uncorrelated  white  noise  of 
covariance  R(K)  .  The  measurement  equation  is  written  as 
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Z(K)  =  h  (X,Kj  +  V(K) 
where  E[U(K)]  =  0 

E[V(i)V(j)T]  =  RS  A  j 

The  initial  conditions  are 
E[X(0) ]  «  X(0)  E[(X(0)  -  XQ)(X(0)  -  XQ)]  =  P0  (4.27) 

E(Wj.Vj^)  =  0  uncorrelated  (4.2S) 

The  Kalman  filter  algorithm  is  to  minimize  the  estima¬ 
tion  error,  i.e.,  error  covariance.  If  X"  denotes  an  esti- 
mate  of  state  vector,  X  ,  and  X  is  the  mean  value  of  the 
state  vector,  the  error  covariance  matrix  is  defined  as 

P  =  E[(X  -  X)  {X  -  X)T]  (symmetric)  (4.29) 

The  development  of  the  extended  Kalman  filter  algorithm  with 
a  few  definitions  is  essentially  more  application  of  the 
expectation  operator  E.  We  will  consider  the  updating 
process  then  the  extrapolating  process.  The  update  equa¬ 
tions  are  used  to  incorprate  the  latest  measurement  in  the 
estimate  and  in  the  covariance.  After  we  obtain  the  updated 
covariance  matrix  we  will  then  consider  the  optimum  choice 
of  Kalman  gain  and  derive  the  equation  for  the  esti¬ 
mation  of  state.  Let  X(K-1)(-)  denote  an  estimate  of  X 


(4 . 25) 

(4.26) 


70 


at  time  K-l  ,  given  measurements  up  to  and  including  ti;:ie 


K-2  ,  P(K-l) (-)  denote  the  error  covariance,  K(K-l) 
denote  Kalman  gain  and  H(k-l)  denote  the  Jaccobean  matrix 
at  time  K-l  .  We  will  seek  an  optimal  observer  in  the 
presence  of  the  state  excitation  noise  and  in  the  presence 
of  the  observation  noise. 


P(K-l)O)  =  E{[I  -  K(K-1)H(K-1)]X(K-1)[X(K-1)TC-)(I-K 
(K-1)H(K-1)T)+V(K-1)TK(K-1)T] 

(4.3 

+  K(K-l)  V(K-l) [  X  (  K  -  1 ) (-)T(I-K(K-1)H(K-1)T) 

+  V( K- 1) TK ( K - 1 ) T ]| 

Rearranging  Equation  (4.30)  using  the  definitions 
Equations  (4.26)  and  (4.29)  and  the  assumption  that  the 
measurement  errors  are  uncorrelated,  then  one  has 


P(K- 1)  (  +  )  =  [I-K(K-1)H(K-1)]P(K-1)(-)[I-K(K-1)H(K-1)T] 


+  K(K-l)  R(K-1)K(K-1)T 


(4.3 


This  equation  updates  the  error  covariance  matrix. 

Then  we  still  need  an  optimal  gain  K.  We  will  choose  to 
minimize  the  diagonal  elements  of  the  covariance  matrix, 
i  .  e  .  , 

J(K-l)  =  E[X(K-1) (+)TX(K-1) (+)]  =  trace  [P(K-1)(+)]  (4.3 


mathematical  model 
system 


Introducing  a  bit  of  matrix  calculus,  one  has 


— -  [trace  A  B  AT]  =  2AB  (4.33) 

aK( 1)'  Cl -K(K-l)HCK-l)  3  =  -H(X-1)T  (4.34) 

Applying  these  formulas  to  Equation  (4.31)  and  solving 
for  K(K-l)  , 

K(K-l)  -  P(K-l)  (-)H(K-1)T[H(K-1)P(K-1)  (-)ii(K-1)T  +  R(K-1)  ]  " 1 

(4.3S) 

This  is  the  Kalman  Gain  Matrix. 

Rearranging  Equation  (4.31) 

P(K-l) (+)  =  [I-K(K-1)H(K-1) ]P(K-1) (-1)  (4.36) 

Then  the  update  state  estimate  at  time  K-l  is  equal  to 
the  extrapolated  state  estimate  at  time  K-2  plus  a  term 
which  weights  the  measurement  residual  via  gain  K(K-l)  . 

X(K-1)(  +  )  =X(K-1)(-)  +  K(K-l)[2(K-l)-h(k-l)  ,^(K-1)  (-))  ] 

(4.37) 

We  need  to  consider  also  the  propagation  of  the  estimate 

A 

X(+)  between  measurements  and  the  propagation  of  the 
covariance  matrix  P(K-1)(+)  between  measurements.  Propa¬ 
gation  of  the  estimate  is  straight  forward  and  involves  only 
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the  application  of  the  state  transition  matrix  and  the 
forcing  terms. 

X'(K)C-)  =  HK-1)X(K-1)  (+)  +  U  ( X  -  1 J  (4.38) 

By  the  definition  (4.29),  solving  for  P(K) (-) 

P(K)(-)  =  i(K-l)P(K-l)  (+)i(K-l)T  +  Q(K-l)  (4.39) 

Figure  4.2  illustrates  the  above  developed  equations  in 
block  diagram  form  which  contains  a  system  model,  measure¬ 
ment  process  and  the  Kalman  filter  to  be  implemented,  and 
Figure  4.3  shows  a  simplified  computer  flow  diagram  of  the 
discrete  Kalman  filter. 

Equations  (4.35)  through  (4.39)  constitute  the  recursive 
formulas  for  implementing  the  extended  Kalman  filter 
algorithm.  The  process  is  initializing  by  providing  values 

A 

x(0)(-)  and  p(0) (-)  .  In  order  for  the  extended  Kalman 

filter  to  be  simulated,  the  next  section  will  describe  the 
initialization  of  the  Kalman  filter. 

E.  INITIALIZATION'  OF  THE  KALMAN  FILTER 

The  one  objective  of  this  work  is  to  test  the  perform¬ 
ance  of  BTT  missile  control  system  having  a  Kalman  filter  as 
an  optimal  estimator  as  a  function  of  different  measurement 
vectors.  All  values  for  initialization  of  former  worker  in 
this  work  (Ref.  14)  will  be  used  for  computer  simulation. 


1.  Initializing  the  Error  Covari  ance  Matrix ,  P 

Let  c  ,  a  and  a  define  the  standard  deviation 
of  relative  position  estimate,  of  relative  velocity  esti¬ 
mate,  and  of  target  acceleration  estimate.  The  initial 
values  of  error  covariance  matrix  p  ,  will  be  chosen  as 
follows.  The  first  six  elements  of  the  initial  state 
vector  are  determined  bv 


X(0)  =  X(0)  +  V(0) 

Where  the  measurement  vector  X  is  the  true  measurement 
vector  X  plus  V  caused  by  measurement  uncertainties . 
The  expected  mean  square  of  initial  covariance  matrix 
components.  With  the  assumption  that  the  initial 
covariance  matrix  is  diagonal  with  the  diagonal  elements 
reflecting  the  uncertainties  associated  with  the  initial 
state  vector  estimates,  the  components  of  the  covariance 


matrix  are 


EtVj2]  =  E  C  V2  2  ]  =  E[V32] 


e[v42]  =  E[VS2]  -  E [ V6 2  3  =  ov2 


The  last  three  elements  of  the  state  vector,  acceleration 
components,  are  set  to  zeros 


Initializing  the  Measurement  Noise  Variance 


Mat  r lx  R 

Measurement  errors  for  the  active  seeker  are 
assumed  to  be  due  to  thermal  noise,  gimbal  angle  pickoff 
error,  environmental  noise,  and  glint.  The  values  of  one 
sigma  errors  in  each  measurement  will  be  used  in  present 
study.  Because  in  the  case  of  angle  measurement,  the  envi¬ 
ronmental  noise  is  proportional  to  the  square  of  range  to 
target  and  glint  errors  the  wander  of  the  apparent  target 
centroid  as  seen  by  the  seeker  as  varing  inversely  with 
range  to  target.  The  errors  caused  by  noise  are  generally 
a  highly  complex  function  of  the  target  geometry,  target 
radar  cross  section,  radar  receiver,  s ignal-to-noise  ratio 
and  etc.  For  simplicity,  in  present  study,  the  constant  one 
sigma  errors  (average  values)  in  each  measured  vector  are 
assumed  to  initialize  the  Kalman  filter.  Referring  to 
Reference  13,  the  one  sigma  values  of  the  average  errors 
for  measurement  of  angle,  angle  rate,  range  and  range  rate 
of  the  present  active  seeker  are  tabulated  below. 


(1) 

angle 

measurement  errors 

:  0.15 

-0.6  deg 

(2) 

angle 

rate  measurement  errors 

:  0.5 

-2.0  deg/sec 

(3) 

range 

measurement  errors 

3  meters 

(4) 

range 

rate  measurement  errors 

6  m/sec 

Before  the  simulation  runs  is  undertaken,  it  will  be 
briefly  discussed  the  manner  in  which  the  Kalman  filter  on 
board  the  missile  is  initialized  at  the  start  of  the 


engagement  based  on  sensor  data  from  the  launch  aircraft. 
The  launch  aircraft's  sensors  are  assumed  capable  of  deter¬ 
mining,  within  some  prescribed  accuracy,  the  relative 
position  and  relative  velocity  of  the  target  with  respect 
to  the  missile  at  the  instant  of  launch.  This  information 
is  sued  to  initialize  the  first  six  elements  of  the  filter' 
state  vector.  However,  the  launch  aircraft  is  assumed  to 
be  able  to  provide  no  information  regarding  target 
acceleration.  Thus  the  three  elements  of  the  initial  state 
vector  are  set  to  zero,  finally,  the  complete  list  of 
nominal  parameters  for  initialization  is  provided  in  Table 
4.1. 

F.  EVALUATION  OF  THE  OPTIMAL  ESTIMATOR  PERFORMANCE 
The  performance  of  the  state  estimator  with  two- 
measurement  vectors  (case  1) ,  with  four-measurement  vectors 
(case  2) ,  and  with  six-measurement  vectors  (case  5) ,  was 
simulated  to  evaluate  the  effect  of  the  possible  implementa 
tion  of  the  more  measurement  sensors  on  the  bank-to-turn 
missile  for  typical  scenario  (tail  chase  engagement  case). 
The  error  covariance  and  the  Kalman  gain  components 
selected,  and  the  states  estimated,  were  computed  as  shown 
in  Figure  4.4  through  Figure  4.30.  The  control  laws  imple¬ 
mented  with  the  estimator  also  were  tested  to  evaluate  the 
performance  of  the  control  system  of  the  BTT  missile  as 
shown  in  Figures  4.51  through  4.35. 


Figures  4.4  to  4.12  represent  the  variations,  of  the 

error  covariances  of  each  state  variable.  Figures  4.1 

through  4.7  show  the  three  covariances  for  the  relative 

positions,  i.e.,  ,  P??  ,  P-_  .  Only  the  covariance 

of  the  two  measurement  case  varies  substantially.  Thus 

one  would  expect  that  the  gain  components  >  G?^  >  G-^ 

would  show  only  variations  in  the  two  measurement  case  as 

are  shown  in  Figures  4.13  through  4.15.  The  trend  noted  in 

the  set  of  graphs  (4.4,  4.5,  4.5)  and  in  the  corresponding 

Kalman  gains  is  also  observed  in  the  other  covariance 

components.  Thus  the  covariance  components  for  the  velocity 

(Figures  4.7,  4.8  and  4.9)  should  larger  variation  is 

reflecting  in  the  corresponding  Kalman  gains  of  4-16,  4-17, 

4-18.  The  covariances  of  the  acceleration  components  show 

a  similar  development  but  there  appears  to  be  less  variation 

in  a  and  acceleration  covariances.  In  the  case  of 

a  the  variation  of  the  Kalman  gains  are  restricted  to  a 

much  narrower  range.  The  error  covariance  and  the  Kalman 

gain  components  in  all  three  cases  eventually  converge  to 

very  small  values.  Figure  4.22  through  Figure  4.30  show  the 

results  of  the  state  estimations  of  the  three  estimators. 

Figures  4.22,  4.23  and  4.24  show  the  variations  of  the 

estimated  relative  position  (R  ,  R  ,  R_) .  Figures  4.25, 

4.26  and  4.27,  the  variations  of  the  estimated  relative 

velocity  (V  ,  V  ,  V  ).  Figures  4.28,  4.29  and  4.30  the 
rx  ry  x  z 
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variations  of  the  estimated  target  acceleration  (a  , 
a.  ,  a.  )  .  The  trends  of  the  curves  reuresenting  the 
estimated  states  are  similar  except  for  Figure  4.25  which 
shows  a  wide  variation  in  the  estimated  V  component  of 
the  two  measurement  case  from  the  true  state  values. 
Generally,  the  estimated  states  of  the  estimator  with 
two-  and  four-measurement  are  close  to  the  true  states, 
in  spite  of  the  w'ide  variation  of  error  covariance  compo¬ 
nents  of  two-measurement  case.  The  estimator  with  six- 
measurement  vectors  generates  the  estimated  states  which 
show  almost  same  characteristics  as  the  estimator  with 
four-measurement  vectors,  but  usually  underestimated  the 
states  in  comparison  with  the  generated  states  of  the 
estimator  with  four-measurement  vectors.  Figure  4.32  and 
Figure  4.54  show  the  curves  of  the  control  commands  com¬ 
puted  using  the  estimated  states.  The  results  of  the  six- 
measurement  case  have  a  wholly  different  form  from  the 
curves  representing  the  computed  control  commands  in  the 
other  cases.  The  normal  acceleration  commands  computed 
using  the  estimted  states  with  the  si x  -measurement  vectors 
initially  have  smaller  values,  but  finally  reach  maximum 
value  as  shown  in  Figure  4.32.  Also  as  shown  in  Figure 
4.54,  the  roll  rate  commands  run  from  the  minimum  limit  to 
the  maximum  limit.  The  large  variations  of  the  both  control 
commands  of  case  3  during  a  short  time  interval  result  in 
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the  worst  miss  distance,  15.4m.  The  control  commands  in  the 
other  cases  behave  almost  like  characteristics  representing 
that  the  control  commands  are  increasing  smoothly  from  the 
zero  to  its  maximum  control  command,  after  two  or  three 
seconds  then  slowly  decreasing  to  converging  values.  The 
relatively  small  variations  of  the  both  control  commands 
during  a  short  time  interval  comparing  with  the  control 
commands  of  case  3  result  in  9.7,  miss  distance  in  case  1 
and  8.7m,  miss  distance  in  case  2.  From  above  results,  the 
estimator  with  four-measurement  vectors  gives  the  best 
result,  the  estimator  with  six-measurement  vectors  the 
worst  result  and  the  estimator  with  two-measurement  vectors 
almost  same  result  as  the  of  four-measurement  case  comparing 
with  the  result  of  the  estimator  with  six-measurement  vectors. 

In  original  measurement  equation,  for  case  1,  the 
missile  has  two-measurement  vectors  which  give  the  informa¬ 
tion  'about  the  relative  angles,  9^  and  ^  ,  was  assumed. 
These  two-measurement  vectors  are  non-linear  functions  of 
three  relative  distance  components,  X^ ,  ,  X,  .  For  case  2, 

the  missile  has  four-measurement  vectors  which  give  the 
information  about  above  two  angles  plus  relative  range  and 
time  rate  change  of  relative  range  which  are  also  non-linear 
functions  of  the  components  of  six  state  vectors,  X,,  X0 , 

X^ ,  and  X^  .  For  case  3,  eventhough  the  missile 

has  six-measurement  sensors  are  assumed,  the  measured 
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vectors  using  the  six-measurement  sensors  contain  only 


state  vectors  which  are  exactly  same  kinds  and  numbers  as 
the  measurable  state  vectors  using  the  four  measurement 
sensors  in  case  2.  As  the  angle  rates,  8D  ,  yn  ,  are 
included  to  the  estimator  as  the  measurement  vectors  in 
case  3,  the  both  angle  rates  are  non-linear  functions  of 

only  six-state  vectors,  X, ,  X-,  X-,  X.,  Xr  and  X,  . 

1  2  a  4  5  6 

The  results  of  the  estimator  with  four-measurement  vectors 
and  with  six-measurement  vectors  show  the  almost  same 
variations  of  the  error  covariance  and  the  Kalman  gain  com¬ 
ponents  of  the  six  measurement  case  as  of  the  four  measure¬ 
ment  case  as  shown  in  Figures  4.4  through  4.21. 

The  Kalman  filter  in  this  work  is  a  first  order  filter. 
Thus  the  components  of  the  Jacobean  matrix  contain  only  the 
first  order  term,  i.e.,  all  higher  order  terms  are  neglected 
As  the  order  of  the  system  is  increased  one  would  expect  a 
more  complex  system  to  response  differently.  The  above 
reasons  affect  to  the  underestimation  of  the  states  of  six 
measurement  case  as  shown  in  Figures  4.22  through  4.50. 

From  the  above  considerations  it  can  be  concluded  that  as 
the  increased  number  of  measurement  vectors  is  implemented 
to  the  extended  Kalman  filter  algorithm,  the  result  of  the 
system  may  not  be  enhanced  up  to  the  expected  degree  due  to 
the  increased  complexity  of  the  system.  The  system  with  the 
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V.  PERFORMANCE  EVALUATION  OF  THE  CONTROL  SYfTE.Y 


This  section  will  discuss  the  performance  of  the  control 
system  implemented  in  a  bank-to-turn  missile  for  each  scen¬ 
ario.  The  maximum  normal  acceleration,  the  maximum  roll 
rate  and  the  time  constant  are  the  key  parameters  which 
reflect  the  missile  capabilities,  the  noises  to  measured 
vectors  are  the  parameters  of  the  environmental  effect,  and 
auxiliary  variables  of  some  importance  are  then  total  engage¬ 
ment  time  and  missile  velocity.  The  mean  miss  distance 
determined  from  the  50  Monti  Carlo  runs  was  used  as  a 
performance  standard  for  the  estimators  of  two-  (case  1), 
four-  (case  2)  and  six-measurement  vectors  (case  3)  on  each 
scenario.  The  simulation  results  are  shown  in  Figures  5.1 
through  5.19.  The  sensitivities  to  variations  of  the  para¬ 
meters  will  be  discussed  below.  Since  the  simulation  results 
show  essentially  the  same  characteristics  for  the  tail -chase 
engagement  (scenario  1)  and  for  the  head-on  engagement 
(scenario  2).  The  performance  on  these  two  scenarios  will 
be  discussed  together  for  each  parameter.  The  side- 
approach  engagement  (scenario  5)  will  be  discussed  sepa¬ 
rately  . 

Figures  5.1  and  S.2  show  the  miss  distance  variation  as 
a  function  of  missile  maximum  normal  acceleration.  For 
scenario  1,  case  1  shows  a  mean  miss  distance  of  9.8m  at 
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above  1 3g '  s ,  case  2.9m  at  above  log's  and  case  5,  15.  2 .z  at 
above  log's.  Each  case  also  shows  that  as  the  maximum  limit 
of  the  normal  acceleration  decreases,  the  mean  miss  distance 
in  case  1  is  less  sensitive  than  in  the  other  two  cases. 

For  scenario  2,  the  result  of  case  1  shows  the  minimum  mean 
miss  distance  is  10m  and  the  result  of  case  2  is  17m  at 
a  =  21g's  respectively,  the  result  of  case  3  is  32m  at 
a  „  =  17g's.  The  trends  of  the  curves  show  that  the  mean 
miss  distances  in  case  1  and  case  2  vary  in  same  manner, 
which  rises  sharply  as  the  maximum  acceleration  decreases. 

On  the  other  hand,  the  mean  miss  distances  in  case  3,  vary 
smoothlv  above  a  .  =  log's  .  One  mav  conclude  from  the 
analysis  as  follows:  for  scenario  1,  the  mean  miss  distance 
can  be  decreased  since  more  information  on  the  state  vectors 
can  be  obtained  from  the  increased  measurement  vectors  in 
estimator  for  low  maximum  acceleration.  In  order  to  get  the 
mean  miss  distance  below  20m  for  all  three  cases,  one  needs 
am„„  0;f  over  llg's.  For  scenario  2,  as  the  measurement  state 
vectors  are  increased  in  estimator,  the  mean  miss  distance 
does  not  decrease,  possibly  due  to  the  effect  of  neglecting 
the  higher  order  terms  in  the  extended  filter.  As  the  mag¬ 
nitude  of  the  values  of  the  missile  and  target  geometry 
components  becomes  larger,  the  effect  of  the  linearization 
becomes  more  pronounced.  amax  =  over  19g's  is  required  to 


to  obtain  the  mean  miss  distance  below  15m  in  case  1  and 
case  2.  In  case  5,  the  mean  miss  distance  dees  no: 
below  20m. 

Figures  S.3  and  5.4  the  variations  of  the  mean  miss 

distances  as  missile  maximum  roll  rate  varies.  In  scenario 

1,  for  case  1  and  case  2,  a  mean  miss  distance  of  below  10m 

can  be  obtained  when  P  is  greater  than  3  rad/sec.  For 

max  ° 

case  3,  one  has  a  mean  miss  distance  of  below  16m  with 
P  greater  than  6  rad/sec.  The  trend  of  the  mean  miss 

distances  is  almost  constant  over  certain  range  in  values 
of  P  „  .  The  variation  of  the  mean  miss  distance  below 
this  range  is  a  tenth  of  the  magnitude  comparing  with  the 
variation  of  the  mean  miss  distances  as  a  function  of 
maximum  acceleration.  It  is  concluded  that  for  scenario  1, 
the  mean  miss  distance  in  case  1  and  case  2  is  not  effected 
by  the  roll  rate  if  the  roll  rate  is  above  3  rad/sec.  In 
case  3,  higher  roll  rate  limit  can  reduce  the  mean  miss 
distance.  For  scenario  2,  the  mean  miss  distances  of  case 
1  and  case  2  are  like  the  trends  for  scenario  1.  However, 
the  mean  miss  distance  of  case  3  has  a  minimum  of  26m  at 
Pmax  *  4  rad/sec,  then  increases  up  to  56m  then  decreases 
slightly.  This  trend  may  be  due  to  the  effect  of  the  roll 
rate  lag,  which  is  caused  by  the  increased  maximum  roll 
rate  and  by  the  underest imat ions  of  the  state  variables  in 
case  3.  The  other  characteristics ,  i.e.,  the  mean  miss 
distance  decreases  with  increased  maximum  roll  rate  limit 
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and  tor  scenario  1,  case  2  shows  less  near  jniss  distant  .• 
tiian  case  i,  but  for  scenario  2,  the  opposite  trend  shivs, 
can  be  explained  with  the  sane  reasons  as  those  of  the  first 
analsv i s . 

The  variation  of  the  mean  miss  distances  with  the  tine 
cosntant  as  a  parameter  is  shown  in  Figures  5.5  and  5.6. 

For  both  scenarios,  one  has  a  flat  plateau  for  low  values  of 
the  time  constant  and  a  sharp  rise  with  increasing  calucs  of 
the  time  constant,  i.e.,  for  scenario  1,  Figure  5.5  shows 
tiie  curves  of  the  mean  miss  distances  slightly  increase  up 
to  the  time  constant  of  0.9  second,  for  scenario  2,  the 
curves  up  to  0.5  second  are  plateau,  then  for  both  scenarios 
the  curves  rise  sharply.  For  scenario  1,  the  mean  miss 
distance  of  case  2  is  less  than  that  of  case  1  up  to  0.S 
second  time  lag,  but  above  0.5  second  lag,  the  results  show 
an  opposite  trend.  Also  over  this  time  lag,  the  curves  of 
case  5  for  both  scenarios  are  not  sensitive  to  the  variation 
of  parameter  up  to  the  time  constant  of  0.7  second  for 
scenario  1,  0.9  second  for  scenario  2.  The  general  trend  of 
the  curves,  i.e.,  as  the  time  constant  becomes  longer,  the 
miss  distance  becomes  larger,  is  due  to  slow  system  response 
From  the  figures  showing  the  results,  the  time  lag  of  less 
than  0.5  second  is  necessary  to  obtain  the  mean  miss 
distance  of  the  required  order  of  magnitude. 

Figures  5.7  and  5.8  show  the  variations  of  the  mean  miss 
distances  as  a  function  of  one  sigma  angle  noise.  For 
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FIGURE  5-6.  MISSILE  TIME  CONST, -NT  VS 
MEAN  MISS  DISTANCE 
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scenario  1,  as  shown  in  Figure  3.7,  the  maximum  miss 
distance  is  obtained  tor  the  six  measurement  case  i.  j  cb  the 
miss  distance  being  the  smallest  for  the  two  measurement 
case.  Above  about  the  angle  error  of  o . 5  degree  the  miss 
distance  is  not  a  strong  function  of  the  angle  error.  For 
small  error  values,  i.e.,  less  than  0.2  degree,  the  curves 
seem  to  be  approaching  each  other.  For  the  scenario  2,  the 
miss  distance  is  a  strong  function  of  the  angle  error.  For 
the  two  scenarios  (Figure  5.7  and  5.S),  the  two  measurement 
case  has  similar  results  up  to  an  error  of  0.4  degree.  All 
cases  show  a  strong  rise  in  the  miss  distance  with  increas¬ 
ing  error.  The  four  measurement  case  seems  pasticularlv 
sens i t i ve . 

The  variation  of  the  miss  distances  as  a  function  of 
missile  velocity  are  shown  in  Figures  5.9  and  5.10.  During 
simulation  the  initial  stand  off  distance  was  increased  to 
keep  the  engagement  time  the  same.  One  important  aspect  of 
Figure  5.9  is  that  represents  an  opposite  trend  of  the 
vatiations  in  the  mean  miss  distance  to  the  general  trend 
of  the  variation  in  the  mean  miss  distance.  Under  the 
missile  velocity  of  700m/s,  the  mean  miss  distance  is  less 
than  10m  for  case  5,  from  15m  to  12. 5n  for  case  2,  from  15.7m 
to  14.5m  for  case  1.  However,  the  curves  for  scenario  2 
show  the  usual  trend,  i.e.,  the  mean  miss  distance  of  case  1 
is  the  smallest  and  that  of  case  3  is  the  largest  among  the 
mean  miss  distances  of  all  three  cases  at  a  value  of  tiie 
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given  parameter.  For  both 

the  curves  slowly  increases  as  the  missile 
increases,  except  tor  one  curve  of  case  1. 
riced  the  foregoing  cons ide rat  ions  as  follows:  \s  : 
values  of  the  missile  and  target  geometrv  components 


":n 


h e c  one 


smaller,  i . e. ,  relatively  less  effect  of  neglecting  the 
higher  order  terms  in  the  first  order  filter  and  less 
contamination  of  the  measured  vectors  to  noises,  the  mean 
miss  distance  can  be  enhanced  with  the  estimator  of  the  six- 
measurement  vectors.  As  the  estimator  algorithm  becomes 
simple,  the  trend  of  the  mean  miss  distances  is  less  sensi¬ 
tive  to  the  target  geometry. 

Figures  5.11  and  5.12  represent  the  effect  on  the  mean 
miss  distance  of  the  total  engagement  time.  For  simulation 
the  initial  stand  off  distance  was  increased  to  keep  the 
constant  missile  velocity.  The  general  behavior  of  the  miss 
distances  of  all  three  cases  slowly  decrease  as  the  engage¬ 
ment  time  increases.  For  case  1,  the  mean  miss  distance 
slowly  decreases  or  is  almost  constant  at  below  10m.  In  the 
other  two  cases,  there  is  more  variation  in  the  miss  distance 
but  the  general  trend  is  still  downward  with  increasing  time 
to  go.  Note  that  as  time  tc  go  increases  for  a  given  T, 
one  has  more  characterist ic  time  interval  available  reflected 

in  T  /t  .  The  mean  miss  distance  therefore  should  be  less 
go 

as  the  total  engagement  time  increases.  This  fact  is  veri¬ 
fied  from  both  figures. 


Figures  5.15  through  5 .  IS  shew  the  variations  e:  ..-.e 
n;ean  miss  distances  of  all  three  cases  for  side -appro  web 
engagement  (scenario  5)  as  a  function  of  each  parameter. 

Tiie  large  miss  distances  resulting  from  this,  scenario  5 
motivated  its  separate  consideration.  General  trends  of  the 
curves  are  similar  to  the  other  scenarios  but  with  higher 
miss  distance.  This  scenario  is  a  strong  test  of  the 
system.  In  this  scenario,  the  control  laws  are  subject  to 
the  error  due  to  the  limitation  of  the  small  angle  assump¬ 
tion  and  the  higher  maximum  acceleration  capability  is 
required  due  to  the  large  variation  of  the  PZC-miss  in 
V-  and  5-direction.  Thus  there  are  some  unexpected  trends. 

Figure  5.15  shows  the  variation  of  the  mean  miss 

distances  for  each  case  as  the  increases.  Although 

the  missile  velocity  is  less  than  that  of  scenarios  1  and  2, 

the  maximum  allowable  acceleration  is  increased.  As  shown 

in  figure,  the  mean  miss  distance  for  case  1  is  11m, 

for  case  2  is  10m,  for  case  5  is  13.5m  only  at  a  x  =  -3g's 

The  figure  shows  the  miss  distance  variations  in  this 

scenario  are  more  sensitive  to  the  a  v.  than  in  the  other 

max 

two  scenarios.  The  mean  miss  distance  of  below  40m  can  be 


obtained  at  a  =  over  19g's,  i.e.,  more  a  ,  is 
max  °  max 

required  than  in  scenarios  1  and  2.  This  trend  of  the  mean 
miss  distances  can  be  explained  with  the  larger  variation  of 


the  state  vector  in  V- direct  ion  which  can  cause  the  .  r .  o 
mg le  U  ?1  in  control  command  equations. 

As  shown  in  Figure  S.14,  the  variations  or  the  tiiss 
distances  of  case  1  and  2  as  a  function  of  p  are  almost 
constant  above  P  ^  =  4  rad/sec.  That  of  case  5,  however, 
continuously  decreases  up  to  P  =6  rad /sec,  then  is 
almost  constant.  For  scenario  3  one  needs  a  maximum,  roll 
rate  of  about  5  or  6  rad/sec  to  have  an  effective  BTT 
missile,  which  is  higher  than  in  scenarios  1  and  2. 

Figure  5.15  shows  the  variations  of  the  mean  miss  distances 
as  a  function  of  time  constant.  Up  to  the  time  constant  of 
0.7  second,  the  mean  miss  distance  is  less  20m.  For  all 
three  cases,  comparing  with  the  results  for  scenarios  1  and 
2,  the  figure  shows  the  curves  of  the  mean  miss  distances 
are  less  sensitive  in  this  scenario.  The  target  accelera¬ 
tion  vector  lies  in  same  direction  as  the  missile  velocity 
as  shown  in  Figure  5.4.  Because  the  components  of  the 
relative  velocity  are  major  components  in  computing  the 
control  commands,  the  more  correct  control  commands  may  be 
computed  due  to  the  small  uncertainty  in  and 

components.  This  results  in  smooth  variation  of  the  mean 
miss  distance  curve.  The  trend  of  the  miss  distances  as 
the  one  sigma  angle  error  increases,  shows  some  different 
characteristics  from  the  trend  for  scenarios  1  and  2  in 
Figure  5.16.  The  mean  miss  distances  for  all  cases 
increase  continuous!}'  as  the  noise  input  increases  up  to  the 


LJislanc 


0.0  0.2  0.4  0.5  0.8  1.0 

Time  Constant  (ssc) 


1.2  1.4 


FIGURE  5-15.  niSSILE  TlilE  COM STRNT  VS 


!EAM  niSS  DISTANCE 


MLssLle  VeLocLty  (m/sscl 


1000.0 


FIGURE  5-17.  MISSILE  VELOCITY  VS 

MEAN  mss  DISTANCE 

144 


SIDE-RPPRCBCH  ENSRSEMEl 

X“u.G5( 1/secl  os“3(mj 

ar.-w  "  .5( dec/sec  j  Pv=y“6.28 !  rat 

9  3  '-s  i  i*  A 

cr.-o-.- .  15 i dec )  c  “23 la's; 

03  w  max  w 

Missile  Velocity  "  BCGlm/secl 

LEGEND 

A  -  2 “MEASUREMENT  VECTORS 
□  -  •» -MEASUREMENT  VECTORS 
O  -  6-MEASUREMENT  VECTORS 


one  sigma  error  of  about  P .  4  degree.  Above  the  one 
error  of  0.-1$  degree  ,  case  1  decreases  then  rises  sh.:r-i’  . 
case  2  continues  decreasing  slowly,  case  5  increases 
continuously.  The  variations  of  the  mean  miss  distances  in 
scenario  1  are  almost  plateau  and  in  scenario  2  are  almost 
constant  up  to  the  one  sigma  angle  error  of  about  0.4  degree 
then  rises  sharply  as  shown  in  figures  5.7  and  5.S.  This 
reason  can  be  deduced  from  the  fact  that  as  the  one  sigma 
angle  error  becomes  large,  the  same  magnitude  of  the  two 
state  vectors  in  different  direction,  i.e.,  relative 
distance  in  X-direction  and  relative  in  Y-direction,  having 
the  same  order  of  magnitude  may  be  contaminated  to  the  noise, 
on  the  other  hand,  the  other  scenarios,  only  the  components 
in  the  X-direction  being  relatively  large  are  contaminated 
to  the  noise.  From  this  fact,  it  can  be  concluded  that  the 
missile  in  a  side-engagement  case  is  more  sensitive  to  the 
angle  error  than  the  other  scenario. 

Figure  5.17  shows  the  variations  of  the  mean  miss 
distances  as  the  missile  velocity  increases  as  a  parameter 
(i.e.,  holding  the  total  engagement  time  of  5  seconds  and 
changing  the  relative  position  of  the  target) .  Up  to  the 
missile  velocity  of  700m/s  the  curves  of  the  mean  niss 
distances  are  plateau,  after  that  velocity  the  curves  rise 
sharply.  Figure  5.8  shows  the  variations  of  the  mean  miss 


distances  as  a  function  of  the  total  engagement  time  (i.e., 
holding  the  constant  velocities  of  the  missile  and  the 
target  and  changing  the  relative  position  of  the  target  in 
scenario) . 

Figure  3.13  shows  the  mean  miss  distances  rise  sharply 
over  the  engagement  time  of  6  seconds.  This  is  somewhat 
unexpected  and  may  be  an  artifact  of  the  implementation 
procedure.  In  other  scenarios  the  mean  miss  distance  as  a 
function  of  the  missile  velocity  increases  slightly  as  the 
missile  velocity  increases  and  that  as  a  function  of  the 
total  engagement  time  decreases  slowly  with  increased 
engagement  time.  From  above  analysis  one  has  a  conclusion 
that  the  computed  control  commands  should  converge  in  the 
final  time  as  the  computed  commands  of  the  two-  four- 
measurement  case  shown  in  Figure  5.32  and  Figure  3.34  to 
obtain  the  low  mean  miss  distance.  However,  it  can  be 
expected  in  this  scenario  that  the  computed  control  commands 
do  not  converge  possibly  due  to  the  violation  of  the  small 
angle  assumption  in  deriving  the  equations  of  the  control 
commands  (I.e.,  in  this  scenario  the  component  is 

large  in  Equation  (5.7),  that  causes  the  large  mean  miss 
distances).  However,  in  the  case  of  relatively  small  magni¬ 
tude  of  the  missile  and  target  geometry  components,  the 
error  caused  by  the  violation  of  small  angle  assumption  is 
relatively  small  comparing  with  the  errors  caused  by  the 
contaminated  of  the  measured  vectors  to  noises. 


limitations  of  the  control  inputs,  end  the  estismz  i 
the  state  vectors  using  first-order  filter.  This  fast  is 
verified  with  the  results  that  the  mean  miss  distance  is 
lira  for  case  1,  10m  for  case  2,  lira  for  case  5  at  the 
missile  velocity  of  6  d  0  m / s  and  the  engagement  time  of  5 
seconds,  as  shown  in  Figures  5.17  and  5. IS. 

The  final  Figure  5.19  shows  the  variations  of  the  mean 
miss  distances  for  case  5  of  all  three  scenarios  as  a  func¬ 
tion  of  the  one  sigma  angle  rate  error.  The  result  repre¬ 
sents  that  the  variations  of  the  mean  miss  distances  due  to 
increased  input  values  of  one  sigma  angle  rate  error  are 
almost  constant  for  scenarios  1  and  2.  For  scenario  3,  the 
mean  miss  distance  varies  at  constant  rate  of  inclination. 
This  slope  of  the  curve  is  relatively  small  comparing  with 
the  slopes  of  the  mean  miss  distance  curves  as  the  other 
parameters  change.  From  this  analysis,  it  may  be  concluded 
that  the  exposure  of  the  measured  vectors  to  the  angle  rate 
error  noise  is  neglegible  in  the  case  of  present  study. 


VI.  CONCLUSIONS 


A  simple  but  effective  biased  guidance  i aw  and  three 
extended  Kalman  filter  equations  have  been  derived  from  the 
optimal  control  theory  for  a  bank- to- turn  missile  with  zero- 
lag  and  with  0.5  second  time  lag  autopilots  for  pitch  accel¬ 
eration  and  roll  state. 

The  equations  of  motion  are  linearized  around  the  present 
orientation  of  the  missile  with  the  assumption  of  small 
future  roll  angles.  During  computer  simulation,  the  optimal 
control  laws  and  estimators  were  tested  separately  invoking 
the  separation  theorem.  In  the  simulation  of  the  system 
three  optimal  estimators  were  used,  i.e.,  estimator 
with  two  measurement  vectors,  with  four-measurement  vectors 
and  with  six-measurement  vectors.  The  system  has  been  eval¬ 
uated  for  three  scenarios,  tail  chase  engagement,  head  on 
engagement  and  side  approach  engagement. 

From  the  analyses  of  the  simulation  results  in  Chapters 
5,  4  and  5,  it  can  be  concluded  as  follows:  The  control  law 
was  successfully  simulated  for  a  hypothetical  bank- to- turn 
missile  with  pitch  acceleration  and  roll  rate  autopilot 
within  the  imposed  limits  on  three  scenarios.  Miss  distances 
with  zero-lag  were  negligible  for  all  three  scenarios 
(below'  0.5m).  For  the  time  lag  case  one  had  miss  distance 
of  5.9m  for  scenario  1,  0.8m  for  scenario  2  and  52.9m,  with 


0.5  second  lag  and  with  a  =  20g's.  Bv  imposin'-  a 
higher  maximum  g  of  23  on  the  scenario  5,  the  miss 
distance  was  reduced  to  below  5m  miss  distance.  One  would 
expect  a  high  g  requirement  for  a  side  engagement  situa¬ 
tion  as  in  the  case  of  scenario  3.  It  is  evident  that  the 
missile’s  tracking  ability  against  the  target  is  very 
sensitive  to  the  missile  and  target  geometry,  the  missile's 
manuevering  capability  is  also  a  very  important  factor. 
Optimal  control  laws  applied  to  the  side  approach  engage¬ 
ment  place  a  severe  strain  on  the  system  from  the  viewpoint 
of  the  small  angle  approximation.  Increasing  the  number  of 
measurements  did  not  result  in  increased  performance  in  the 
sense  of  smaller  miss  distances  for  the  missile  except  in 
isolated  cases.  This  may  be  due  to  the  increased  com¬ 
plexity  of  the  system  as  the  number  of  measurements  was 
augmented.  Further,  in  this  study  only  first  order  extended 
filters  were  implemented.  It  is  possible  that  better 
performance  could  be  obtained  at  the  expense  of  increased 
computational  load,  by  using  a  second  order  extended  Kalman 
filter. 


APPENDIX  A:  PROGRAM  LISTING 


This  appendix  provides  listing  of  the  computer  program 
used  in  the  present  study.  Only  one  program  for  the  six 
measurement  case  is  provided.  Except  for  a  few  difference 
in  the  subroutine  filter  the  programs  for  the  two  other 
cases  are  identical  to  the  program  included  here. 
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TT  <2CC)  ,AAC  (2  Jj)  ,iPC  i  2C0)  ,  6AAK  (200)  ,  ?J1  (200) 
A  AH  (200)  ,  ETZ  I2..2) 

-Oil  |2C0)  .0  21  (20C)  ,G3  1  (2  30)  ,G4  1  (200)  ,  05  1  (200) 
Cl  ( 2 CO )  ,02  1 100)  ,Ui  (ICC)  ,04  (2CC)  ,05  (200) 

0  6  2 CO)  ,C7  (200)  ,Ue  (2  C0|  ,09  ( 20 C) 

21  (2  C  C)  ,22  1200)  ,23  (200)  ,24  LoO)  ,25  (200) 

£6  (2C0)  .17  ( iO  0 )  .Zi  i  2C  0)  .23  (20C) 

PI  (2CC)  ,P2  (230)  ,P3  UCO)  ,?4  (200)  ,  P5  (200) 
r1.. :  i  n  i  i  o  4  .  ?  r  ■**  *  -a  u  .  ~  -  r  \ 


o. 


ISEU1  DAI  A 

CT1=.0S 
h  C  A3  £=  1 
N  CAX  =  1 10 
KE  EINT  =  1 
LE2XNI=5 
mc:=  i 
IHlOaO. 
ESIO=Q. 

?  HI  0=0  . 

V £0=10  00  . 
I.'1C=0. 

rr.o=o. 

zkc=o. 

110=2500. 

IIC=0. 

ZTC=-100. 

y ixc=5  c  o . 

VI  10=0. 
VI2Q=0. 

E  1=.  00  3254 
E  2  =  1 . 2  5 
iZ'.(0=0. 
AIV0=4 .*5.6 
ATZJ=0. 

Cl  A  =  .  0  5 
CIV=. 0  5 
CIZ  =  . 05 
SAI=5. *9.6 
0=1. /2U. 
f  2  1  =  .  0  026 
R22=.0C37 
£23=. 0  02  £ 
H2«=.uC67 
£25=3. 

326=6 . 

SE  1=.0026 
SH2=.0C87 
3  p  3  =  .  0  0  2  6 
S3«=.OCd7 
SE5=3. 

S36=6. 

SECS=3. 

SVEL=6. 

St CC=9.*S.8 
IS?IS«.5 
T  ACC=. 5 
ASAX=2C.»S.8 
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uuu 


non 


Q  (5,6)  ='2  2 3 
t  i*.?)  =C-!-l 

115111(6,5) 

.31X3 (6 . c) 

EC  13*.  1=1,3 

■*rx:s  jb,7)  j  (i.i)  ,;  (1,2)  ,c  (i  .3) (i,  2)  ,c  (1,5) ,  iii,o) ,;  (i,  7) 

*«;!;)  hen.') 

CC.STIN  uc 
?C?.3AI  (•  1*1 

FC5S.M  (//SaX,  'CUCISiL  0  3X7311*//) 

FCKMSI  (55  1-..3) 

rCeilAT  <//4**X,  *  C  .'.A7SXX  2  IIH  PS2UEC  30X55,  iP3U  =  ’  ,  75 . 3//) 

AEC  c3  5U  X  C  NCI35  TO  J  (lAISIX 

I11=(«'?5E**2)  *  |EXT>*-)/4. 

?12=<3?SC**2)  »  <EIX»'*3)/2. 
f  25=  (Ji?SU««2>  *  (E  17**2) 

C  I M)  -Cl J,1)  *511 
Q  U,2)  *1  I*,-)  *5  11 
C  IJ.3)  -1  U,3)  *5  11 

0  I?,?)  «■ 1  *522 

(.  (3,0)  =  („  (5,3)  *1<t 

V  (c,o)  =u  (6, a)  +  522 
Q  ll.?)  =Q  (!»?)  *512 

(,  (4,3)  =>2  (4,3)  *514 
v  U,°)  *1  (3,o)  *5  12 

y  (1,1)  =0  (1,1)  *512 
C  (5,2)  ='2  (5.2)  *51* 

Q  (5,3)  =Q  (c , J )  *514 

KP.IT2J6, 5)  jESC 
EC  133  1=1,3 

1iEI73  j 6 , 7)  2  (X  ,1)  (X.2)  ,i(X,3)  ,;  IX,  a)  ,;  (1,5)  ,  3(1,6)  ,2  (1,7) 

CO  150  1*1,9 

cc  mo  j= i , ) 

cc sxxacs 
saxi: (6,5) 

SXJBT  CA55  LOCE 
DO  950  111=1 , KCA35 
INI71A  LI Z5 
H  =  C 

K=((EaiST-1 
L  =  Lr?.IBl-1 
ISTC?=  0 
01=017 

IF  (3CA  £2.50.1)  GC  XC  160 
XJEEI'.CT  =  1  CCoO  J 
IFSX.'JT=  1  CCCOO 
I f 1CT=  0 
K  =  C 
1=0 

CCSIIM02 

I  (1)  =0. 

I  (2)  =THXC 

x  (3  ?es;c 


I  (2)  =THIC 

?(4=»si5( 

x  (5)  =x:iu 
r  (6 )  =  x  no 
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nnn  nnn  nnn  n  n  nfir,  nnn 


non  n  non  o  o  n  n  non 


122=031:1  irHi)  *csi:.  iz.iT)  *Diia  (isi)  ■ 

a22«3si:«  iehii  *:coi  i : 

J31OC03  uHIl  *0SIX  UH7)  *OCOS  U3t» 

232-sccs  i * Hi)  »:su  i ::-.i)  *ssia  usu 

133=CCC3  ifHI)  *£C0i  (TEI) 

VSX=A1  1  *  VC 
VX7  =  A  1  2*  VC 
V  £2=  11  3*  VS 

CC22UTS  TfJS  CCXEC.'iFNIS  OF  5X113  ' 

aXsTT-XX 
8  3  =  31-  i.'l 

sa*zr-ia 

VFX=7?.C-V»X 
V?.3  =  VT3-  \  S  t 

vf:=vto-vc  z 


«C„05  ( J : 
‘  ox  i:i  ( ?  5 : 


n:(  =  ATX0*33X3  (-CXX*T1 
Aai»i;*C*C2i?  <-CI'i«T) 
AT2-AI  20  ♦  CL.(2  (*C.T'*='j 

X  2  ( 1 )  =  3  X 
X I X  2:  =8  3 

XI  (3)  =  SZ 
X2  (4)  =  VS  ;t 
X 2  (5)  =  V  3  3 
2 <  <3>  =  V 5  2 


X-  IT  =  AT .< 

:<4  (3i  =ir! 

X2  |9)  =  AI2 

XGC  =  -  (  FXCVFX*-?  3  3*7  32)  / <7=X*= 

FCSJ  3JASCS33  -SI  7  3CIC.’. 


f:<3-*2) 


6£IS*3 

SGI  33  =  (3  X«73X  ♦  £S*'»n  3*53’  VS 


7RXJ-  (C3X.V  (THT3)  *DSI  S  (8  5X3)  *SY* 
35  IS)  »Si*733)  )/  (S*S«GCCS  (THIS)  ) 


ISIIIALZ2E  ST  A  1 3  7EC703 

IF  Cl.  11.  1)  3C  TC  750 

CALL  L'.iC  i  C  1 1  5  3  5  3  ,  >  •  ,  c  ,  1  ,  1  ) 

X  (1)*3X*::-L£  !  .  ■«  (  !)  )  *  f  C  5 
X  1 1)  «S  3*  Or  ,3  I*  *  US  )  *3r03 
X  (3)  «ai*  r  Li (••  (J))  *3?05 
X  (4)»/sX»;: Li  («•  (4)  !  ’S73L 
X  (5)  =733*C3L3  I  »'•  I  5,1  j  *37  3  i 
X(?i)=7SO*LLL2(.'<mo)J  *37  3L 
X  I")  =u  . 

X  IS  I  =  3. 

X  IS)  =0. 

CCMIHUE 

£11=0003  (THTS)  »CtCS  (JilS) 

C  »2»OCOi  li:u  3)  •Oil:.  (  »3I3J 
013**0311 (THIS) 

X£CS=D  1  1*X  (  1)  *£  K»X  (2)  *0  1  >X  l  5) 


non  non  n  n  nnn 


>3113,' 6,0  0  1,F  n,i'2,?33,?44,F55,--5c,F'7f  .-SS,^ 

FCS2AT  (/  1CF13.3) 

FCSM AT ( tOf IS. 3) 

CC  VIIN CE 

iFji-ia.iE3i.Na)  >jc  ac  3 sc 
5BITS(6,2C)  a 

DC  i«  I=1.J  . . , 

W5i:S(6f  j  C)  (If  1 )  i<inl  ^{1 1  1)  <I#i)  /  j  .i  i  •'  i  -  *  •)  »  jit  i.<  ( t#  ■ ) 

*G  A  Zli  J  X  f T  —  ' 

?cSAAi'(//5FjC,  "1AIJ  2.ATSEX  it  =  •  ,F5 . 2  ,  1  )  1  /) 

icsxiri'px, .->:*!  2.=) 

V»  Vrf  »»  *T  a  .»  ‘Jm 

CHECK  FCF  TSSSIJAIICX 

I  r  ( SO  T  .  3  E  -  C .  3  .  C  j  -  GGC  ■  L.E  .  3  •  3 )  GC  TC  AGO 
IF  iIGG.G1.GG)  GC  GC  4..0 
C I *150 

is:cp=  1 
CCMDIUS 

IHIEGHAIE  IHAGECTC3Y  CHE  STEF  AHEAD 

CAii  ::;teg  (•*,  =  > 

IF  (13X0?. 31.0)  GC  1C  700 
IF  (K.uc.  if  AX)  CC  10  SCO 

EXECUTE  XAiKAS  FIITEF  (UFCAXE  STATE  YEC1QE) 

CAII  FILTxE(X,F,E,G,GAi:;) 

GC  TO  70 C 

TEEMIMiL  ESIM  A  S.  C  F  ICT 
CCSTIHUE 

CZ  =  Y}lij -?B) 

CM*DS0i!T  lCX*»2»Ci«''2*02*»4) 

X  2IS5  (II I )  =  c :< 

Y  2 133  (III)  s0  Y 


CALI 


EX.CY.EZ.DM.TTI 

1  =  1  ,?  iC.J.SX,  'OX* 
X  ,  ■  213  S  =  '  ,?  10.  3  ,5 

', PIC  .3.51,  '3*' 

=  '  ,F10.3,5X 

X,  'CASE-  '  ,17) 

GC  TC  950 

El, 5,  J.1,1, 'TI2E 

(SEC)  ',10,'SX 

(MST333) ' 

C1,s,f  .»,0,'TI2E 

(SEC)  ',10,'SX 

(MET20S)  ' 

>.  ’»  'TIME 

(SEC)  •  ,  1C,'  21 

iil.!.,2,1,0, 'TIME 

1 .  ,5. ,5.) 

3 3  1  ,1 , 1  ,  'TIME 

(SEC)  •  ,10,'SY 

(3ETE5S)  • 

(SEC)  ',10, 'EE 

(METEES)  ' 

,C3,a,2,I,3,'TIS£ 
i  ....  G>) 

(SEC)  '  ,  1  0  ,  '  EG 

(METES S)  ' 

,34,3,1, 1,  1  ,  'TIME 

(SEC)  ',  1C,  'V?.X 

(M  ETE3S/SEC) 

,U4  , <  ,  1,  0,  'TIME 

(SEC)  '  ,  1C,'  VKX 

(METE  33/SEC) 

.  D  .  ,5 . 1 

.  15,5,1,1,  1,  'TIME 

(SEC)  '  ,  1C,'  73 Y 

(METE  iS/  ESC) 

:T?i.N,’3;\,o,  'time 

(SEC)  '  ,  10,  '  VSY 

{MSTE3S/S  EC) 

)  .  ,  3  .  ,  5  .) 

~  AJ  1  1  1  l"**s 

(SEC)  ',10,  '732 

(MEXESS/SEC) 

159 


55 

6C 


,  16,  J. , C.  ,0.  , 

CALL  ?tCTG(r7 

1c,  0. 


CALL  PLOT  G  ("1,27 ,6,1  ,1,  1  , 

*  1  ^  I  0  •  *0*  f«;« 

£All  PLCIGJTI^LO,!*,*,!,  j, 

#  1  7  ,  J  •  / 


CALL  PL  _ 

^  19^0*  /  0 •  / 'j »  / 
C  A 


CALL  ? LCTo  f  XI  ^  Lc  f  li ,  ,  1  r  J  # 

t  a  t  LU? Lot g  (i : ",  1 5 ' *i ,  i , 

,  19,0, 

CALX.  ?  LC 10  I  I  *  ,  c 

/Cm  /3m  ,  3  m  (J  ■  /2m  1 


■C.,0. , 


,0. , J.  ,  0 .  , ’J 


\ 


■  C  •  ■  0  •  ,L*,v«,3«,3a) 

Lai  PLCZG(?:.ttc.:i.:,i,i 
.•U03.1AL  ACC2l25»-t:CU  (.1/ 
OIL  PLCTG/TT  i  ,  J 


i ,  1. 

*»  1 

J, 

♦TIM  £ 

(SEC) 

’ , 1C,' 

78  Z  (3STE3S/SEC)  ' 

I-1- 

1  , 

♦TIME 

(SEC) 

'.10,' 

ATX  (.TEIESS/CEC**2) 

%'S, 

3. 

'TIME 

(SEC) 

' ,  10, ' 

4IX  (M  IT  Z  3S/SEC®  *  2) 

=  .) 

1,1. 

1  , 

'TIME 

(3SC) 

•.10,  ' 

AT  t  { M  ST  L  HS/  3  ZC  *  *  l) 

?•) 

*  ,  1  , 

3  .) 

0, 

'TIME 

(SEC) 

■,1C,* 

ATY  (:iETIr;j/3£C*^2) 

1,1. 

1, 

I  ’li’u  ; 

(SEC) 

’,10,’ 

ATT  (MSTI 35/ JEC*  *2) 

s  ,  1 , 

3.) 

0, 

♦TIME 

(SEC) 

• , 10.  1 

ATT  (:iST2iS/5HC**2) 

1  .1. 

1 , 

»*T  -#  T 

(SEC) 

MO, 

♦S^2T  (FI  1)  '  ,  3 

1.1. 

1 . 

'TIME 

(SEC) 

MC, 

' S QHT ( ? 22 )  '  .9 

1.1. 

1, 

'TIME 

(SEC) 

•,1C« 

'  S  2?.  T  ( ?  3  3 )  ’,3 

•  *  4  t 

1, 

*T  IM  E 

(SEC) 

■ ,  1C. 

<S«SI(?94)  ',<3 

1,1, 

1 . 

'TIME 

(SEC) 

'  ,  13. 

•SQHT  (?55)  '  ,  9 

1  .1. 

1 , 

•  'T^  V  R 

(2  SC) 

',10, 

‘5  2 2. 1  (=66)  '  ,  3 

1  .1, 

1 . 

'TIME 

IS  EC) 

MO, 

'$  521  (177)  •  ,  9 

1.1. 

1, 

•TIME 

(SEC) 

MC, 

*3«a:(?fl8)  '  ,  3 

1.1, 

1, 

'TIME 

(SEC) 

MC, 

•555I(S9»)  ’,2 

95  C 


,  1  MJ:U  nl  ACCEibaA.Lo:.  i  ../ 

CALL  ?L03G<II  .BBC  , ’,1  ,  1 
,  1  FCLL  SATE  (aAC/SEC)  ',  19 
CALL  ?[/0ZG(7Z,;Bc,::,£.1.3 
SAXE  (EAt/SECj,  19 
CALL  ?  LCTG  (TL  ,  2  P3S  .  S  ,  1.  1 , 

,  ‘EA.HX  A:;GL2  iCZG)  f  ,  16, 0. 

call  plctg<i:,tg2,:;,  1,1,1 
,  '2I!15  LC  '3Z  |  EEC)  •  ,  Id,  5. 

CALL  ?L0?G(T?,IG1,a,;,l. 0 
.'THIS  TC  90  (SEC)  lb,  j. 

CALL  ?  LC1G ( IT  ,S1  I  ,9,1.1  ,  1 
,'G(2.2)  ■,b,0.,0.,0.,J.,5 
CALL  PLOTG  (IT  ,G2  1  ,:.,1,1  ,  1 
,'Gj5,a)  l,b,0«,C,,C,,v*,S 
CALL  PLCTG  (XT  ,  Gj  1,S,  1,1,  1 

fci!Jd&LkT6j**;«rJ;Ci;i 

till*  *  VczW-'/tv'.r.rt'.i 

*  '  G  (S  f  1)  '  /  o  f  C «  i  u  «  ■  0  •  f  i/>  r  3 

Cali  plcho.  ,cl,9=9) 

CCMTIli  CE 

C0KPUIE  KISS  CIEIASCE  STATISTICS 

17  (NCASE-EO.  1)  GC  TC  975 

CALL  SEAT  (X.TISS,  .‘  CAS  E  ,X3  A a,  X  PI’S) 

CALL  STAT  ( i  HISS,  1«  CASE,  X  0  AE ,  Yaf.S) 

CALL  STAT  (3.'*  ISS  ,  1.  C  AS  i  .  3  3  A  P ,  3  ?  .»£ ) 

CALL  STAT  (OXI  SS  ,  *C  AS E ,03  AS,  0 B  JS) 

9  SITE (6,51 
KBITS  (6,5  5) 

KBITS  (6 , 6C)  XEAf.XEXS 
W  HITE  (  6 , 7  C)  TcAS.JBKS 
K  BITE  (  6  ,  cC)  3EAB,:3!'S 
WRITE (  6 , 9Cl  3EA5,tBES 

?C  B  MAT  (/  /  13X  ,  1  X  la  S  CIST  .VICE  STATISTICS'/) 
7C3KAT  {/5X,'  X-CC:iEC3ESI'  ,5X,  >  ?.EA:.=  ‘  ,  rb.  2,5i, 


/3m/0m/2m/Z 


.  'TIKE  (SEC)  •  ,  10 
;.-.C  “  •;)  '  ,  ,0  ,  j .  ,  9  . 

,  '  T I  *4  '  ( S — C  )  *  ,  1 0 

03 3  AL  _  ACC  Tib  A  A IIC 1.  iX./3EC**2)  l,30,C«,C.,G.,0«,5,,5.) 

*  ‘  71 1!  ^SSC)  ' ,  10  . 

■  ,  V  •  ,  J  •  /  3  .  ,5.  ,5.) 


■  THE  (SZC)  ’,10 
," . , C,,u,,0  •  ,5,  ,  5  . ) 
1,  ’-mir.i  (SEC)  ’Jo 
0.  ,  C.,G.f5*,5«) 

#T IMS  (S£C)  » e  10 

'l{EE  (stc)  1 1  10 

0.  ,  C.f  J.,5  ..D.) 
"*1.VE  (5Z.C )  *  ,  10 

’  T 1  i  E  (SEC)  ’,13 

*  T 1 1  E 

<hh 

<liU 

.  3.  > 


(SEC)  •  ,  10 
(SEC)  • ,  10 
(SEC)  ’  ,  10 


1113=  *  ,  F3. 2) 


[ail 


SOEHC'OTIMH  ?I  ETca  (X  ,1,3,2  ,o) 

IXPLEL’ET  5£.\l*E  ;.x-H,C-i5 

eisessicn  a r> ,si  ,  at  <■/,»)  ,x<0)  ,sr)) ,?  |9t3»  ,iO,9)  ,ax (  <i.  , 

DI2S:iS!CS  ASM  |5, S)  ,2X{  3)  ,??(-;  ,-•)  ,n  (  t ,  j )  ,  a  I  i  ■)  ,  o(  ,  ..  :i  -  i  ' ,  -  I 

C 2 2c. NS  ICS  a? S'!  (6  ,c)  ,3  io  ,  c  )  ,  I  10  ,o  )  ,  3  {  ,  ,a  )  ,  ari  :  ->)  ,  i o)  ,  {-;  j 

CI2£j3iCl  EX  (3)  ,3  rJ  ,  3J  ,  30  A)  ,0ci.(  3  ,  3)  ,  2 1  It  7  (o  ,-,j  ,  <w:,  (■>) 

CCWOS/CSc/  A  ,  AT  ,£:L  ?. 

SXIHAP0LAT3  STATE  VECT03 

CATE  V2'Jlf?(A  .  J,3, 9,  1,-3, 3,AX,  9,11) 

DO  SO  1=1,3 
XX  II)  =  AX  li)  -S  (1) 

cc:;tx:ius 

EX22APOLA7C  3  C  li  A  A  2  A  A  C  c 

c.ul  vaci?r{r  ,ai,s,3,9,9,9.?ai,j,i2) 

CALL  VSCLEF  (A  ,?*T,S,S,3  ,9  ,3  ,  A  EAT,?,!  3) 

DC  15  J  1  =  1, S 
EC  100  J  =  1,9 

P  ?  IX, J)  =AFAI  1 T  ,«,}  ♦*  (X,w) 

CCSTaLu  c 
CCNTIMUS 

Fcaa  a  a:ic  :ii  tatsices 

X1  =  XX<1) 

X2=XX  (  2) 

X  3  =  X  X  (  3  J 
X4=XX  (4) 

X5=XX  5 
X6=XX  (6) 

D  *DSQS?  |X1**2*X2«*2*:<3**2) 

D1=DS 0»T  <1)**2*X2*»2) 

0  2=  X 1 »  X  4  ♦  X  2  *  X  6  *  X  3  *  X  o 

E1=X4»Xo 

E2=X5+  X6 

E3»X4-X5 

E4=D1*Q1*2.*D*C 

D3*.X1*X3*  (2-  »E*C*C1«t1«i  1-  (X  1*X1*E1  »X2*X2*E2>  «  (3.  «0  1  ‘‘J  HD*E)  ) 
C4=X2*  X3*  (2.  *E*D*C  1*C1*E2-(X  1**1  *£1*  X2*X2»E2)  *  (J.  *3  1*0  1*0*D)  ) 
ES=  (XI  *X1*E1*  J2»  X2*E 2)  »  (D»D-X5*X3) 

26=X2*S3*  f0*E*C1*C1*X1»X1*E4) 

07=X1*S3*(-0*D*D1*C1*X2*X2«E4) 

CC  250  1=1,6 
EC  200  J*1,9 

a  (i.Ji  =o. 

CCXTISUS 
CCNTIU  UE 


H  (  ,1' 
B  <1,2 
a  ii,3 
H  (2,  I 
H  ic, 2 

a  <2, 3) 
a  (x,4) 

H  (c ,  S | 
H  <c,6) 
Hli'l) 

a  <3,2! 

H  (4,1) 
a  (4,^) 
a  <4,3) 

118:3 

a  <s,  i» 

K  15,2) 


=  X  1**3/ <C1*C*C) 

=  X  2  *X  3/  1  C1*X*C) 

=-C1/  D  *  C) 

=03/  <  0**5)  »  1C1*=3)  ) 

=  04/  (i  D»*5l  »  (C  1*«  J)  ) 

=  0 S/I  0**5)  »C1) 

=  X  1  *X  1  *  X  3/  lC*E*D*01) 

=  x 2 *.<2 * x 3/  (:*:  *:*o  ij 
=  0  1*X3/1C*C*E) 

=-X2/(C1*C1) 

=  X1/<0  1  *  C  1 ) 

=  Ct/(<C*»3)»  (C1*=4)  ) 

=  C7/  I1c**j(  * ( C  1  *  *  4 )  ) 

=<i1*X2*>3«:3)/  (<0*«3)  *  |E1««2)  ) 

=  X  I/O 
=  X  2/  0 
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